DTIC 

iELECTEJ 
JANO  139111 


MICROMECHANICAL  modeling  of  the 
FIBER/MATRIX  INTERFACE  REGION  IN  A 
METAL  MATRIX  COMPOSITE 
THESIS 

David  D  Robertson 
Captain,  USAF 

AFIT/GAE/EN  Y/90D-22 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 


AIR  FORG 


iTITUTE  OF  TECHNOLOGY 


U| 


AFIT/GAE/ENY/90D-22 


) 


MICROMECHANICAL  MODELING  OF  THE 
FIBER/MATRIX  INTERFACE  REGION  IN  A 
METAL  MATRIX  COMPOSITE 
THESIS 

David  D.  Robertson 
Captain,  USAF 

AFIT/GAE/ENY/90D-22 


Approved  for  public  release;  distribution  unlimited 


AFIT/GAE/ENY/90D-22 


MICROMECHANICAL  MODELING  OF  THE  FIBER/MATRIX 
INTERFACE  REGION  IN  A  METAL  MATRIX  COMPOSITE 


THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 
In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Aeronautical  Engineering 


David  D.  Robertson 
Captain,  USAF 


December  1990 


Approved  for  public  release;  distribution  unlimited 


Acknowledgements 


I  express  my  deepest  thanks  to  Dr.  Shankar  Mall  for  his  calm  patience 
and  expertise  in  guiding  me  through  this  effort.  His  many  insightful  sug¬ 
gestions  and  enlightenments  were  invaluable  and  provided  the  necessary 
catalyst  to  keep  me  on  the  right  path  in  spite  of  my  unwary  attempts  to  the 
contrary.  In  addition,  I  express  my  appreciation  to  Mr.  Ted  Fecke  and  the 
Turbine  Engine  Division  of  V/ right  Research  and  Development  Center  for 
their  interest  in  sponsoring  this  tilers. 

Furthermore,  I  would  like  to  thank  Captain  Howard  Gans  for  his 
extensive  instruction  in  MSC/NASTRAN  and  SDRC-IDEAS  and  for  the 
use  of  an  excellent  computer  system.  Also,  I  am  indebted  to  Dr.  Anthony 
Palazotto  to  whom  I  owe  my  initial  understanding  of  the  finite  element 
method.  Likewise,  I  express  my  gratitude  to  all  the  users  of  lincoln  who 
exhibited  great  patience  as  I  devoured  memory  and  processing  time,  and 
thank  you  Randy  Aldridge,  a  friend  of  mine  and  employee  of  Apple 
Computers  Inc.,  for  helping  me  obtain  a  first  class  home  computer  at  a 
fraction  of  the  cost.  No  doubt,  I  would  still  be  trying  to  write  this  if  I 
owned  an  IBM. 

Finally,  I  reserve  my  most  grateful  applause  for  my  wife,  Diane.  Her 
understanding  and  hard  work  have  been  an  inspiration  and  made  this  epi¬ 
sode  in  our  life  a  pleasant  one,  and  may  God  bless  two  small  children  with  ’  - - 

W 

a  Daddy  who  was  forced  to  refuse  many  invitations  to  play  with  them.  □ 


ii 


a 


i  By _ 

:  Distribution/ 


Availability  Codes 


Dist 


Af-iil  and/or 
Special 


Table  of  Contents 


Page 

Acknowledgements . ii 

List  of  Figures .  v 

List  of  Tables . ix 

List  of  Symbols . x 

Abstract .  xiii 

I  Overview .  1 

1 . 1  Introduction . 1 

1.2  Background .  3 

1.3  Proposal . 4 

II  Perspective . 5 

2.1  Motivation . 5 

2.2  Micromechanical  Methods . 6 

2.3  Imperfect  Bonding .  10 

III  Methods  of  Analysis .  19 

3.1  Hardware  and  Software .  19 

3.2  General  Information  on  Models  Used  in  this  Study .  23 

3.3  Finite  Element  Models  for  Strong  and  Weak  Interfaces.  .  .  25 

3.4  Models  for  Incorporating  an  Interphase  Zone .  30 

3.5  Combined  Interphase  Zone  and  Weak  Bond  Models .  39 

IV  Results  and  Discussion . 44 

4.1  Strong  and  Weak  Interfacial  Bond  Results . 44 


iii 


«.V» 1  • 1  'u.Ui. ~ .1  m«. .^.>.1 , j.j. ■ ... !.,. ■ . • .  ^  j.-u1 


Page 


4.1.1  Results  from  Perfectly  Strong  Interface .  45 

4.1.2  Results  from  Perfectly  Weak  Interface .  50 

4.2  Results  from  Models  Possessing  an  Interphase  Zone .  66 


|  4.2.1  One-Dimensional  Results  With  Interphase  Zone.  .  .  66 

4.2.2  Finite  Element  Solutions  With  Interphase  Zone.  ...  75 

4.3  Results  of  Combined  Interphase  Zone /Weak  Bond . 86 

Models 

4.3.1  One-Dimensional  Results  of  Combined  Interphase.  .  88 
Zone  /  Weak  Bond  Model 

4.3.2  Finite  Element  Results  of  Combined  Interphase  ...  94 
Zone  /  Weak  Bond  Model 


V  Conclusions .  101 

VI  Recommendations  and  Suggestions .  108 

Bibliography .  110 

Appendix  A:  Sample  MSC/NASTRAN  Data  Deck .  114 

Appendix  B:  Mathematica  Program  for  a  Simplified  One- .  117 

Dimensional  Micromechanics  Model  With  an 
Interphase  Zone 

Appendix  C:  Mathematica  Program  for  a  Simplified  One- .  121 

Dimensional  Micromechanics  Model  With  a 
Combined  Interphase  Zone  /  Weak  Bond 

'vita . 128 


iv 


List  of  Figures 


Figure  Page 

1  Nimmer’s  Results  for  Strong  vs.  Weak  Interface,  Plane  ....  15 
Stress  Finite  Element  Solution 

2  Nimmer’s  Results  of  Analysis  vs.  Experiment,  3-D .  16 

Generalized  Plane  Strain  Finite  Element  Solution 

3  Basic  Cell  Used  in  Analysis . 24 

4  35%  Fiber  Volume  Fraction  Square  Array  Finite  Element  ...  26 

Grids  for  Beth  Strong  and  Weak  Interfacial  Bonds 

5  32%  Fiber  Volume  Fraction  Rectangular  Array(Aspect . 29 


Ratio=1.2)  Finite  Element  Grids  for  Weak  Interfacial  Bond 

6  Simplified  One-Dimensional  Model  Incorporating  an . 31 

Interphase  Zone 

7  Coarse  and  Fine  Grids  for  Interphase  Zone  Size  5%  of . 37 

Fiber  Radius 

8  Finite  Element  Grids  for  Interphase  Zone  Sizes  of  7.5%  ....  38 
and  10%  of  Fiber  Radius 

9  Simplified  One-Dimensional  Model  With  Combined . 40 

Interphase  Zone  /  Weak  Bond 

10  Interphase  Zone  Location  Scheme  for  Combined  Finite  ....  43 

Element  Model 

1 1  Strong  Interface  Plane  Stress  and  Plane  Strain  Finite . 46 

Element  Model  Residual  Stress  After  Cooldown 

12  Strong  Interface  Generalized  Plane  Strain  Finite  Element.  ,  .  47 

Model(3-D  Model)  Residual  Stress  After  Cooldown 

13  Strong  Interface  Transverse  Response  for  Plane  Stress .  48 

Plane  Strain,  and  Generalized  Plane  Strain  Finite  Element 
Solutions 

14  Schematic  of  Out-of-Plane  Constraints  Showing  Difficulty.  .  49 

of  a  Plane  Strain  Solution  in  the  Presence  of  Large  Thermal 
Strains 

15  Weak  Interface  Transverse  Response  of  RVE#1  Assuming  .  .  51 

Linear  Elastic  Matrix  and  Using  a  Plane  Stress  Finite 
Element  Solution 


v 


I 

r 

> 


r 


r 

[ 

> 


H 


I 


i 

> 

l 

S 


Figure  Page 

16  Weak  Interface  Deformed  Geometry  of  RVE#1  at  140  MPa  .  .  53 
and  210  MPa 

17  Weak  Interface  Transverse  Response  of  RVE#1  With  and  ...  54 
Without  Friction  at  the  Interface  and  Assuming  a  Linear 
Elastic  Matrix 

18  Weak  Interface  Transverse  Response  of  RVE#1  Assuming  .  .  55 
Elastic-Plastic  Matrix  With  a  Yield  Stress  of  800  MPa  and 

No  Interfacial  Friction  as  Compared  to  Nimmer’s  Results 


19  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#1 . 56 

With  Elastic-Plastic  Matrix  and  Frictionless  Interface  at 
Cooldown  and  Onset  of  Fiber/Matrix  Separation 

20  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#1 . 57 

With  Elastic-Plastic  Matrix  and  Frictionless  Interface  at 

210  MPa  and  260  MPa  Transverse  Load 


21  Weak  Interface  Transverse  Response  for  RVE#2  Using  ....  60 

Plane  Stress  and  Generalized  Plane  Strain  Solutions  With 
Comparison  to  Nimmer’s  Results 

22  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#2  at.  .  .  61 

Maximum  Load  for  Both  Plane  Stress  and  Generalized 

Plane  Strain 

23  Weak  Interface  Transverse  Response  for  RVE#2  Plane . 63 

Stress  Solution  With  and  Without  Fiber/Matrix  Slip  at 

the  Interface 

24  Weak  Interface  Transverse  Response  for  RVE#2  General-.  .  .  64 
ized  Plane  Strain  Solution  With  and  Without  Fiber/Matrix 

Slip  and  Comparison  to  Nimmer’s  Experiment 

25  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#2  at.  .  .  65 
Maximum  Load  Using  Generalized  Plane  Strain  Solution 

and  Preventing  Fiber/Matrix  Slip 


26  One-Dimensional  Solutions  With  Interphase  Zone  for .  68 

Interphase  Yield  Stresses  from  75  MPa  to  200  MPa  as 
Compared  to  Nimmer’s  Strong  and  Wes.k  Bond  Model 

27  One-Dimensional  Solutions  of  Interphase  Zone  Strain . 71 

Hardening  Parameter  Variation 

28  One-Dimensional  Solutions  of  Interphase  Zone  Size .  72 

Variation 

29  One-Dimensional  Solutions  of  Interphase  Young’s . 73 

Modulus  Variation 


vi 


Figure 

Page 

30 

One-Dimensional  Solutions  of  Interphase  Zone  Thermal.  .  . 
Expansion  Coefficient  Variation 

.  74 

31 

Finite  Element  Plane  Stress  Solutions  of  Interphase  Zone.  . 
Yield  Stress  Variation  as  Compareu  to  Perfectly  Strong 
and  Perfectly  Weak  Bonds 

.  76 

32 

Von  Mises  Stress  Contours  After  Cooldown  Using  Plane  .  .  , 
Stress  Solution  for  Strong  Interface  and  for  an  Interphase 
Zone  Yield  Stress  of  175  MPa 

.  78 

33 

Von  Mises  Stress  Contours  After  Cooldown  Using  Plane.  .  . 
Stress  Solution  for  Interphase  Yield  Stresses  of  125  MPa 
and  75  MPa 

.  79 

34 

Finite  Element  Plane  Stress  Solutions  of  Interphase  Zone.  . 
Strain  Hardening  Parameter  Variation 

.  80 

35 

Finite  Element  Plane  Stress  Solutions  of  Interphase  Zone  .  . 
Size  Variation 

.  81 

36 

Finite  Element  Generalized  Plane  Strain  Solutions  of . 

Interphase  Zone  Yield  Stress  Variation  as  Compared  to 
Perfectly  Strong  and  Perfectly  Weak  Bonds 

.  82 

37 

Finite  Element  Plane  Stress  Solutions  of  Conventional  and. 
Fine  Grids  Using  Baseline  Parameter  Set  of  Table  5 

.  84 

38 

Comparison  of  One-Dimensional  and  Finite  Element . 

Transverse  Response  for  Baseline  Parameter  Set  of  Table  4 

.  85 

39 

Comparison  of  Nonlinear  Characteristics  Between  the  .... 
One-Dimensional  and  Finite  Element  Solutions 

.  87 

40 

One-Dimensional  Solutions  of  Combined  Interphase  Zone  . 

/  Weak  Bond  for  Weak  Bond  Ratios  from  0.0  to  1.0 

.  89 

41 

One-Dimensional  Solutions  of  Combined  Interphase  Zone  . 

/  Weak  Bond  With  Loading  Sequence  Explanations 

.  90 

42 

One-Dimensional  Solutions  of  Combined  Interphase  Zone  . 

/  Weak  Bond  for  Various  Weak  Bond  Ratios 

.  91 

43 

One-Dimensional  Solutions  of  Combined  Interphase  Zone  . 

/  Weak  Bond  With  a  Weak  Bond  Ratio  of  0  50  for  Various 
Interphase  Zone  Yield  Stresses 

.  92 

44 

One-Dimensional  Solutions  of  Combined  Interphase  Zone  . 

.  93 

/  Weak  Bond  With  a  Weak  Bond  Ratio  of  0.90  for  Various 
Interphase  Zone  Yield  Stresses 

vii 


Figure 

Page 

45 

One-Dimensional  Solutions  of  Combined  Interphase  Zone  . 
/  Weak  Bond  for  Various  Strain  Hardening  Parameters 

.  94 

46 

Finite  Element  Plane  Stress  Solutions  of  Combined  Inter-.  . 
phase  Zone  /  Weak  Bond  Models  for  Various  Weak  Bond 
Ratios  -  High  Range 

.  95 

47 

Finite  Element  Plane  Stress  Solutions  of  Combined  Inter-.  . 
phase  Zone  /  Weak  Bond  Models  for  Various  Weak  Bond 
Ratios  -  Low  Range 

.  96 

48 

Finite  Element  Plane  Stress  Solutions  of  Combined  Inter-.  . 
phase  Zone  /  Weak  Bond  With  a  Weak  Bond  Ratio  of  87% 
for  Various  Interphase  Zone  Yield  Stresses 

,  .  98 

49 

Finite  Element  Generalized  Plane  Strain  Solution  of . 

Combined  Interphase  Zone  /  Weak  Bond  at  7%  Weak 

Bond  Ratio  as  Compared  to  Plane  Stress  Solution 

.  99 

50 

Finite  Element  Generalized  Plane  Strain  Solution  of . 

100 

Combined  Interphase  Zone  /  Weak  Bond  at  87%  Weak 
Bond  Ratio  as  Compared  to  Plane  Stress  Solution 


List  of  Tables 


Table  Page 

1  Material  Properties  Used  in  me  Finite  Element  Solutions  ...  15 
of  Figure  1 

2  Material  Properties  Used  for  Strong  Bond  Models  and . 45 

Weak  Bond  RVE#1  Models 

3  Material  Properties  Used  for  Weak  Bond  RVE#2  Models  ...  59 

4  Baseline  Parameters  for  One-Dimensional  Model  With . 67 

Interphase 

5  Baseline  Parameters  for  Finite  Element  Models  With . 75 

Interphase 

6  Summary  of  Interface  Models  and  Their  Characteristics.  .  .  105 


ix 


List  of  Symbols 


a  Fiber  Dimension  of  Simplified  One-Dimensional  Models 

a  Thermal  Coefficient  of  Expansion 

(Xf  Thermal  Coefficient  of  Expansion  for  Fiber 

(Xi  Thermal  Coefficient  of  Expansion  for  Interphase 

(Xm  Thermal  Coefficient  of  Expansion  for  Matrix 

b  Interphase  Dimension  of  Simplified  One-Dimensional 

Models 

P  Symmetric  Angle  Defining  the  Weak  Bond  Portion  of  the 

Interface  for  the  Finite  Element  Combined  Interphase 
Zone  /  Weak  Bond  Model 

2 p/90  Weak  Bone  Ratio(P  in  degrees)  for  Finite  Element  Com¬ 

bined  Interphase  Zone  /  Weak  Bond  Model 

c  Matrix  Dimension  of  Simplified  One-Dimensional  Models 

d  Weak  Bond  Dimension  in  Combined  Interphase  Zone  / 

Weak  Bond  Simplified  One-Dimesional  Model 

d/a  Weak  Bond  Ratio  for  One-Dimensional  Combined  Inter¬ 

phase  Zone  /  Weak  Bond  Model 

E  Young’s  Modulus 

Ep  Young’s  Modulus  of  Fiber 

Ei  Young’s  Modulus  of  Interphase 

Em  Young’s  Modulus  of  Matrix 

Z  Total  Strain 

£F.  Eli ,  £12  Total  Strain  in  Regions  F,  Ij ,  and  12  in  Simplified  One- 

Dimensional  Models 


i>  £M2’ 

£m3 


Total  Strain  in  Regions  Mi,  M2,  and  M3  in  Simplified 
One-Dimensional  Models 


x 


eEi1 »  £Ei2 

Elastic  Strain  in  Regions  Ii  and  I2  in  Simplified  One-Di¬ 
mensional  Models 

8Pl1>  £Pl2 

Plastic  Strain  in  Regions  Ii  and  I2  in  Simplified  One-Di¬ 
mensional  Models 

ep 

Load  Error  Fraction 

ew 

Energy  Error  Fraction 

f 

Coefficient  of  Friction 

H' 

Strain  Hardening  Parameter 

V 

Poisson’s  Ratio 

RVE 

Representative  Volume  Element 

RVE#1 

Representative  Volume  Element  for  35%  Fiber  Volume 
Fraction  and  Square  Fiber  Array(Aspect  Ratio=1.0) 

RVE#2 

Representative  Volume  Element  for  32%  Fiber  Volume 
Fraction  and  Rectangular  Fiber  Array(Aspect  Ratio=1.2) 

0 

Stress 

OF,  CJip  0 12 

Stress  in  Regions  F,  1 1 ,  and  I2  in  Simplified  One-Dimen¬ 
sional  Models 

OMj,  Om2, 
OM3 

Stress  in  Regions  Mi,  M2,  and  M3  in  Simplified  One-Di¬ 
mensional  Models 

Oa 

Applied  Transverse  Load 

Oys 

Yield  Stress 

oysj 

Initial  Interphase  Yield  Stress 

ays^ 

Present  Yield  Stress  in  Region  Ii  Accounting  for  Strain 
Hardening  in  Simplified  One-Dimensional  Models 

Oysi2 

Present  Yield  Stress  in  Region  I2  Accounting  for  Strain 
Hardening  in  Simplified  One-Dimensional  Models 

0  x,  0  y,  0  z 

Normal  Stresses  Along  Axes 

XI 


0  VON  MISES 

Von  Mises  Stress 

T 

Temperature  After  Cooldown 

Tp 

Processing  Temperature 

Txy.  Tyz, 

?zx 

Shear  Stresses  Parralle!  to  Faces  Defined  by  Axes  x,  y, 
and  z 

u 

Displacement 

UF>  Uli,  UI2 

Displacement  of  Right  Hand  Faces  of  Regions  F,  1 1 ,  and 
12  in  Simplified  One-Dimensional  Models 

UMi-  UM2- 
UM3 

Displacement  of  Right  Hand  Faces  of  Regions  Mi ,  M2, 
and  M3  in  Simplified  One-Dimensional  Models 

UG 

Displacement  of  Left  Hand  Face  of  Region  Mi  in  Simpli¬ 
fied  One-Dimensional  Model  with  Combined  Interphase  / 
Weak  Bond 

Vf 

Fiber  Volume  Fraction 

AFIT/GAE/ENY/90D-22 


Abstract 

i 

in  order  to  develop  reasonable  analytical  models  and  to  better  under¬ 
stand  how  the  integrity  of  the  fiber/matrix  interface  in  high  temperature 
titanium  based  metal  matrix  composites  effects  the  overall  transverse  nor¬ 
mal  characteristics,  a  parametric  study  of  various  fiber/matrix  interface 
conditions  was  performed.  Analysis  methods  included  both  a  simplified 
one-dimensional  mathematical  model  and  a  finite  element  analysis  using 
MSC/NASTRAN.  First,  both  strong  and  weak  interfacial  bonds  were 
examined.  Subsequently,  to  provide  for  partial  bonding,  a  model  that 
used  an  elastic-plastic  interphase  zone  between  the  fiber  and  matrix  was 

examined. -^Finally,  a  combined  model  was  investigated  which  assumed 

Up 

that  a  portion  oFthe  interface  behaved  as  a  *  /eak  bond  and  the  remainder 
as  if  an  elastic-plastic  interphase  zone  existed  between  the  constituents. 

Results  from  a  weak  fiber/matrix  bond  indicate  that  interfacial  separa¬ 
tion  and  slip  have  a  major  impact  on  the  response,  resulting  in  substantial 
matrix  plastic  flow.  At  the  point  of  separation  a  significant  nonlinearity 
in  the  transverse  response  curve  appears,  so  that  the  modulus  after  separa¬ 
tion  is  considerably  weakened.  Slip  along  the  interface  was  also  found  to 
significantly  effect  the  transverse  response,  and  increasing  the  interfacial 
friction  arrested  the  slip  only  before  separation  had  occurred.  Therefore, 
after  separation,  no  capability  existed  for  varying  the  transverse  response. 
An  additional  weak  bond  model  was  created  that  allowed  only  fiber/matrix 
separation(no  slip  was  allowed).  The  additional  constraint  of  no  slip  was 
found  to  greatly  stiffen  the  response. 

xiii 


Since  a  weak  bond  model  offered  no  capability  for  varying  the  level 
of  bonding,  the  elastic-plastic  interphase  zone  model  was  incorporated  to 
determine  if  partial  bonding  could  be  obtained.  It  was  found  that  by  vary¬ 
ing  the  plastic  properties(i.e.  yield  stress  and  strain  hardening)  any  level 
of  interfacial  bonding  was  possible,  from  perfectly  strong  to  perfectly 
weak.  The  elastic-plastic  interphase  zone  model  also  allows  for  energy 
dissipation.  However,  for  low  interphase  yield  stresses,  a  relaxation  of 
the  thermal  residual  stresses  occurs.  The  combined  interphase  zone  / 
weak  bond  model  arrested  this  relaxation  of  thermal  stresses,  but  did  not 
offer  as  much  flexibility  for  varying  the  level  of  bonding. 


MICROMECHANICAL  MODELING  OF  THE  FIBER/MATRIX 
INTERFACE  REGION  IN  A  METAL  MATRIX  COMPOSITE 


I  -  Overview 


1.1  Introduction 


For  over  two  decades  micromechanical  analysis  of  unidirectional  com¬ 
posites  has  proven  to  be  a  useful  tool  in  predicting  the  overall  properties 
of  composite  lamina  and  in  understanding  how  the  fiber/matrix  interaction 
effects  these  overall  properties.  Correct  modeling  of  the  fiber/matrir  in¬ 
terface  i«  a  crucial  factor  in  a  micromechanical  analysis.  Traditionally, 
for  polymer  or  epoxy  based  composites  a  perfectly  strong  bond  is  assumed 
at  the  interface.  A  vast  amount  of  research  has  been  performed  on  these 
types  of  composites,  and  as  a  result  a  large  body  of  knowledge  exists  to 
confirm  the  strong  bond  assumption.  Unfortunately,  such  materials  are 
only  useful  for  low  temperature  applications.  Therefore,  all-metallic 
structures  have  remained  the  obvious  choice  for  structures  subjected  to 
high  temperatures. 

It  has  long  been  a  desire  to  incorporate  the  benefits  of  composite  ma¬ 
terials  in  high  temperature  applications.  This  desire  has  intensified  as  a 
result  of  programs  such  as  the  National  Aerospace  Plane  (NASP)  and  the 


Integrated  High  Performance  Turbine  Engine  Technology  (IHPTET)  Ini¬ 


tiative.  A  key  element  in  the  success  of  these  programs  is  the  effective 
implementation  of  high  temperature  materials  with  vastly  improved  capa- 
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bilities.  This  technology  thrust  for  advanced  high  temperature  materials 
has  fueled  an  interest  in  other  classes  of  composites  than  the  traditional 
polymer  based  composites.  Among  these  are  fiber  reinforced  metal  matrix 
composites  which  are  considerably  attractive  from  the  manufacturing 
technology  and  material  performance  aspects.  In  particular,  systems 
using  titanium  or  titanium  aluminides  as  the  matrix  and  silicon  carbide  for 
the  fiber  are  appealing  due  to  their  ability  to  maintain  high  strength  and 
stiffness  at  temperatures  in  excess  of  800°C  without  the  severe  weight 
penalties  associated  with  most  existent  high  temperature  metals. 

High  temperature  metal  matrix  composites  present  a  new  challenge  in 
micromechanical  analysis.  The  metal  matrix  typically  possesses  a  much 
higher  thermal  coefficient  of  expansion  than  the  fiber,  and  in  addition, 
processing  temperatures  are  necessarily  much  greater  than  those  used  for 
conventional  thermosetting  composites.  Such  a  combination  results  in 
relatively  h  rge  thermal  residual  stresses  in  the  composite  after  cooldown 
to  room  temperature.  Also,  at  the  fiber/matrix  interface  a  reaction  zone  is 
present  where  the  fiber  has  penetrated  into  and  reacted  with  the  matrix  re¬ 
sulting  in  a  third  phase  material  which  is  likely  to  possess  distinct  proper¬ 
ties  from  either  the  fiber  or  matrix.  Further,  it  has  become  apparent  from 
experiment  that  the  fiber/matrix  contact  cannot  be  modeled  as  a  strong 
bond[l].  Such  conditions  necessarily  make  transverse  micromechanical 
modeling  of  high  temperature  metal  matrix  composites  a  complex  task. 

The  extent  to  which  these  unique  conditions  may  effect  the  overall 
material  properties  are  not  clear  at  present.  Moreover,  in  order  to  develop 
reliable  micromechanical  models  for  titanium  based  metal  matrix  compos¬ 
ites  thorough  investigations  must  be  performed  that  will  ascertain  which 
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characteristics  need  to  b..  included  and  how  they  should  be  modeled. 

1.2  Background 

Performing  a  micromechanical  analysis  of  transverse  normal  loading 
in  a  unidirectional  composite  while  incorporating  both  a  weak  fiber/ma¬ 
trix  bond  and  thermal  residual  stresses  requires  a  nonlinear  material 
analysis.  Residual  compressive  stresses  exist  at  the  interface  after 
cooldown,  and  must  be  overcome  before  the  matrix  will  separate  from  the 
fiber.  Thus,  before  separation  the  composite  behaves  as  if  a  strong  bond 
exists  at  the  interface  while  afterward  it  demonstrates  weak  bond  charac¬ 
teristics.  This  was  experimentally  observed  by  Johnson  et.  al.[l].  Also, 
Nimmer[2;  3]  analytically  examined  this  effect  using  both  a  one-dimension¬ 
al  square  fiber  model  and  a  nonlinear  finite  element  analysis. 

Substantial  energy  dissipation  occurs  during  a  loading  and  unloading 
sequence.  This  energy  loss  is  most  likely  due  to  mechanisms  at  the 
fiber/matrix  interface.  The  only  means  of  energy  dissipation  in  Nimmer’s 
model  were  transverse  friction  along  the  fiber/matrix  contact  and  plastic 
flow  of  the  matrix.  Since  much  of  the  energy  loss  in  Nimmer’s  model  oc¬ 
curred  in  the  form  of  matrix  plastic  flow,  it  exhibited  significant  perma¬ 
nent  strain  after  unloading  which  has  noi  been  observed  experimentally. 

In  addition,  the  model  demonstrated  a  weaker  secondary  stiffness  that  oc¬ 
curs  after  fiber/matrix  separation  than  demonstrated  in  experiment. 

Thus,  the  problem  suggests  that  there  are  additional  constraints  at  the 
interface  that  cannot  be  modeled  using  a  perfectly  weak  bond  assumption 
alone.  If  accurate  load  histories  of  the  composite  are  to  be  reproduced, 
these  additional  mechanisms  must  be  incorporated. 
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1.3 


Proposal 


To  develop  accurate  models,  parametric  studies  must  first  be  per¬ 
formed  to  see  how  various  interface  conditions  effect  the  overall  compos¬ 
ite  characteristics.  The  present  study  will  examine  transverse  normal 
loading  using  a  variety  of  interface  conditions  and  seek  to  model  the  ener¬ 
gy  dissipation  mechanisms  using  a  third  zone.  First,  an  examination  simi¬ 
lar  to  Nimmer’s  [2; 3]  will  be  undertaken  analyzing  the  effect  of  a  strong(no 
separation  or  slip)  versus  weak(both  separation  and  slip  allowed)  interface 
on  the  matrix  stress  and  overall  characteristics.  Comparisons  will  be 
made  with  Nimmer’s  results  for  validation.  Second,  the  drawbacks  of  a 
completely  weak  interface  model  will  be  examined  through  adding  an  ad¬ 
ditional  constraint  that  prevents  slip  at  the  fiber/matrix  interface.  Third¬ 
ly,  a  reaction  zone  at  the  interface  will  be  added  to  model  the  three  phase 
system.  The  zone  will  possess  different  material  properties  than  the  fiber 
or  matrix  in  an  attempt  to  model  the  energy  dissipation  mechanisms 
present  near  the  fiber/matrix  interface.  Finally,  a  portion  of  the  reaction 
zone  will  be  removed  to  model  a  partially  bonded  interface. 

The  problem  at  hand  is  a  highly  nonlinear  one.  Therefore,  even 
though  a  simplified  one  dimensional  model  will  also  be  used,  the  mode  of 
analysis  will  mainly  be  the  finite  element  method.  The  finite  element  pro¬ 
gram,  MSC/NASTRAN,  is  a  highly  versatile  code  that  possesses  nonlinear 
capability  and  a  gap  element  essential  for  modeling  contact  type  problems 
such  as  a  weak  fiber/matrix  bond. 

This  investigation  will  demonstrate  possible  models  and  discuss  their 
credibility.  Such  studies  will  advance  the  knowledge  base  of  high  temper¬ 
ature  metal  matrix  composites,  and  lead  to  their  ultimate  application. 
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II  -  Perspective 


2.1  Motivation 

Before  a  new  material  may  be  incorporated  into  structural  designs  it  is 
essential  to  be  able  to  reliably  predict  its  useful  life  under  a  given  load¬ 
ing.  The  mechanisms  involved  in  fatigue  failure  operate  on  the  micro¬ 
scopic  level.  Therefore,  before  such  predictions  are  possible  a  thorough 
understanding  of  the  micromechanics  involved  is  necessary. 

As  mentioned  in  the  introduction  the  recent  technological  push  for  im¬ 
proved  high  temperature  materials  has  resulted  in  new  classes  of  compos¬ 
ites  being  brought  to  the  forefront  of  composite  research.  Titanium  based 
metal  matrix  composites  have  the  potential  to  meet  this  new  challenge, 
but  they  possess  unique  characteristics  which  must  first  be  understood. 
Among  these  :  e  the  imperfect  fiber/matrix  bond  and  the  significant  ther¬ 
mal  residual  stresses  after  processing. 

Many  questions  arise  when  studying  such  phenomena.  For  instance, 
what  level  of  debonding  occurs  at  the  interface  and  in  what  directions,  or 
how  much  relaxation  of  the  thermal  stresses  occurs  during  both  cooldown 
and  afterward?  Most  importantly,  how  can  these  effects  be  modeled  and 
which  ones  are  significant?  Before  such  questions  may  be  answered,  more 
research  must  be  performed,  both  analytically  and  experimentally.  The 
goal  being  to  develop  proper  micromechanical  models  which  will  lead  to  a 
clearer  understanding  of  the  behavior  of  high  temperature  metal  matrix 
composites. 
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2.2 


Micromechanical  Methods 


The  work  presented  in  this  paper  will  analyze  transverse  normal  load¬ 
ing  using  a  simplified  one-dimensional  model  and  a  numerical  finite  ele¬ 
ment  approach.  However,  many  micromechanical  theories  have  been  pro¬ 
posed  over  the  years,  and  it  is  helpful  to  present  a  brief  history  for  the 
reader.  Also,  before  attempting  a  micromechanical  analysis  it  is  benefi¬ 
cial  to  possess  an  overall  view  on  both  closed  form  and  numerical  meth¬ 
ods  that  have  been  used  in  the  past. 

Chamis  and  Se:ideckyj[4]  render  a  comprehensive  critique  on  various 
techniques  in  micromechanics.  They  divide  the  approaches  into  many 
classes:  netting  analyses+,  mechanics  of  materials  approaches,  self-con¬ 
sistent  models,  variational  techniques  using  energy  bounding  principles, 
exact  solutions,  statistical  approaches,  discrete  element  methods,  semiem- 
pirical  approaches,  and  microstructure  theories.  Many  of  these  approach¬ 
es  are  summarized  by  Jones  [5:85-146].  Except  for  the  first  two,  all  ap¬ 
proaches  use  at  least  some  of  the  elements  of  the  theory  of  elasticity. 

Only  a  few  will  be  briefly  mentioned  here. 

The  self-consistent  model[6]  is  constructed  of  a  fiber  and  matrix  mod¬ 
eled  as  concentric  cylinders  embedded  in  the  composite.  In  this  method 
the  isotropic  plane  strain  elasticity  solution  in  cylindrical  coordinates  is 
solved  for  a  transverse  applied  stress  in  the  surrounding  material.  Such 
an  approximation  is  advantageous  because  it  results  in  a  closed  form  elas¬ 
ticity  theory  solution  allowing  stresses  at  the  fiber/matrix  contact  to  be 
analyzed.  Therefore,  it  has  been  widely  used.  However,  it  completely  ig- 

+  Netting  analysis  assumes  the  fibers  provide  all  the  longitudinal  stiffness  and  the  ma¬ 
trix  provides  all  the  transverse,  shear,  and  Poisson’s  effect. 


6 


nores  the  effects  of  adjacent  fibers(assumes  no  packing  array). 

Other  elasticity  solutions  in  micromechanics  include  closed  form  se¬ 
ries  solutions  such  as  developed  by  Sendeckyj[7]  for  longitudinal  shear 
loading.  Nevertheless,  the  application  of  such  solutions  can  be  tedious 
since  it  involves  the  truncation  of  an  infinite  series. 

Also,  a  very  useful  set  of  equations  for  overall  stiffness  were  devel¬ 
oped  by  Halpin  and  Tsai [8].  Their  result  represents  a  good  approximation 
to  the  self-consistent  model  in  the  transverse  direction  without  the  com¬ 
plexity.  Overall  stiffness  properties  may  be  obtained  quickly  and  easily 
from  relatively  simple  equations.  In  the  Halpin-Tsai  equations  a  single 
parameter  is  used  to  incorporate  fiber  geometry,  packing  geometry,  and 
loading  conditions.  The  value  of  this  parameter  must  be  assessed  by  re¬ 
lating  the  equations  to  exact  solutions.  Also,  no  fiber  or  matrix  stresses 
are  determined,  only  overall  stiffness  properties. 

Except  for  the  self-consistent  model,  micromechanical  analysis  gener¬ 
ally  requires  assuming  a  regular  array  of  fibers.  After  assuming  a  fiber 
array,  symmetry  is  used  to  produce  a  representative  volume  element  con¬ 
sisting  of  a  single  fiber  and  its  surrounding  matrix.  In  doing  this  the  ana¬ 
lyst  has  taken  a  completely  random  distribution  of  fibers  and  assumed  it 
to  be  ordered.  The  influence  of  random  filament  packing  on  transverse 
stiffness  was  examined  by  Adams  and  Tsai[9].  They  studied  two  types  of 
random  arrays:  square  and  hexagonal.  The  arrays  were  not  entirely  ran¬ 
dom  since  they  still  contained  repeating  units  of  several  fibers.  Regard¬ 
less.  results  indicated  that  a  random  hexagonal  array  analysis  agrees  more 
closely  with  experiment  than  a  square  random  array  analysis.  This  con¬ 
firms  the  notion  that  a  hexagonal  array  seems  more  physically  realistic. 
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On  the  other  hand,  if  an  analysis  that  is  not  random  is  used,  the  square 
fiber  array  tends  to  agree  more  closely  with  experiment.  In  addition,  a 
square  array  possesses  greater  flexibility  by  allowing  for  two  independent 
dimensions  in  defining  a  unit  cell  v/hile  the  hexagonal  array  contains  only 
one  dimension.  Therefore,  since  most  analysis  assumes  a  repeating  regu¬ 
lar  array  of  fibers,  a  square  array  is  more  widely  used. 

Much  work  has  been  done  using  numerical  approaches  to  obtain  dis¬ 
crete  approximations  to  the  elasticity  equations,  and  since  the  study  pre¬ 
sented  in  this  paper  will  use  a  numerical  finite  element  technique,  a 
strong  review  of  these  methods  is  in  order.  Early  on  in  1966  Foye[lO]  used 
a  finite  element  approach  with  a  regular  square  array  of  fibers  to  analyze 
the  transverse  properties  of  a  unidirectional  composite.  Subsequently, 
Adams  and  Doner[ll;  12]  analyzed  longitudinal  shear  as  well  as  transverse 
normal  loading  of  a  unidirectional  composite  using  a  finite  difference 
technique.  Stress  concentration  and  stiffness  were  analyzed  as  functions 
of  fiber  volume  fraction,  constituent  stiffness  ratio,  and  fiber  shape. 
However,  the  method  allowed  for  only  linear  elastic  analysis. 

Later  Adams[i3]  used  a  nonlinear  material  analysis  to  analyze  trans¬ 
verse  normal  loading.  A  plane  strain  finite  element  scheme  was  incorpo¬ 
rated  to  obtain  the  stress  field  throughout  the  matrix,  but  the  greatly  in¬ 
creased  computational  requirements  resulted  in  only  a  limited  amount  of 
data  being  produced.  However,  substantial  matrix  yield  was  shown  to  in¬ 
duce  only  a  mild  nonlinearity  in  the  composite. 

Others  sought  to  incorporate  inelastic  micromechanics  through  numer¬ 
ical  techniques.  Foye[l4]  examined  longitudinal  shear  as  well  as  trans¬ 
verse  normal  loading  for  carbon/epoxy,  boron/epoxy,  and  boron/aluminum 
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using  a  nonlinear  finite  element  program.  Again,  the  computational  re¬ 
quirements  limited  the  analysis  to  an  extremely  coarse  grid  and  sparse 
data.  Nevertheless,  it  demonstrated  the  capability  of  inelastic  microme¬ 
chanics  using  the  finite  element  technique.  Later,  Crane  and  Adams  [15] 
performed  a  more  extensive  inelastic  analysis  of  longitudinal  shear  load¬ 
ing  of  a  unidirectional  composite 

Seeking  to  advance  the  capability  of  computational  micromechanics, 
in  1973,  Adams [16]  presented  a  rudimentary  crack  propagation  scheme  to 
model  failure  in  the  matrix  around  the  fiber  duri  .g  transverse  loading. 
Once  an  element  reached  its  critical  octahedral  shear  strain  for  failure,  its 
stiffness  was  reduced  to  zero.  Of  course,  such  a  scheme  is  not  entirely 
accurate  since  it  produces  relatively  large  gaps  in  the  material  after  fail¬ 
ure  of  a  matrix  element.  Also,  modeling  crack  propagation  in  such  a  dis¬ 
crete  fashion  resulted  in  a  highly  irregular  transverse  stress-strain  re¬ 
sponse.  In  spite  of  such  approximations,  the  two  points  on  the  stress- 
strain  curve  associated  with  directly  before  an  element  failure  and  directly 
after  tended  to  bound  the  experimental  data,  but  not  in  a  way  that  could  be 
easily  predicted.  In  addition,  these  earlier  nonlinear  analysis  studies  were 
only  capable  of  monotonically  increasing  loads. 

Continued  work  in  the  micromechanics  field  has  resulted  in  many  var¬ 
ied  prediction  techniques  available  for  today’s  engineer.  The  underlying 
assumption  in  most  of  these  is  a  strong  fiber/matrix  bond.  Also,  only  a 
few  methods  are  able  to  account  for  the  thermal  residual  stresses  due  to 
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Such  limitations  produce  no  difficulty  for  epoxy  based  com¬ 


posites  that  possess  apparent  strong  bonds  and  are  processed  at  relatively 
low  temperatures.  However,  as  previously  mentioned,  such  assumptions 
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are  grossly  in  error  for  high  temperature  titanium-based  metal  m  rix 
composites.  Therefore,  if  rei  able  micromecharical  models  are  to  be  de¬ 
veloped  for  such  composites,  an  imperfect  bond  must  be  incorporated. 

2.3  Imperfect  Bonding 

An  interphase  zone  between  the  constituents  is  known  to  exist  in  all 
composites.  This  zone  comprises  imperfections  due  to  absorption  and  re¬ 
action  of  the  two  constituents,  voids,  imperfect  adhesion,  and  microc¬ 
racks.  Since  the  load  transfer  between  fiber  and  matrix  depends  on  the 
strength  of  this  interphase  zone,  the  thermomechanical  properties  of  the 
composite  are  strongly  affected  by  its  integrity.  However,  due  to  high 
temperature  processing  and  the  constituents  involved,  the  effect  of  the  in¬ 
terphase  zone  in  high  temperature  metal  matrix  composites  is  amplified. 
Therefore,  even  though  an  accurate  representation  of  the  zone  is  very 
complex,  incorporating  its  effect  into  a  micromechanics  model  is  neces¬ 
sary  for  such  composites. 

Various  attempts  have  been  made  to  model  an  imperfect  fiber/matrix 
bond.  Several  researchers  have  used  linear  elastic  methods  while  more  re¬ 
cently  some  inelastic  models  have  been  proposed.  In  the  elastic  realm, 
three  studies[17-19]  employed  the  finite  element  analysis  with  an  interfa¬ 
cial  layer  between  the  constituents.  The  elastic  properties  of  the  layer 
were  varied  to  model  various  levels  of  debonding.  Later,  Lene  and 
Leguillonf20]  allowed  tangential  slip  at  the  inte  face  using  a  finite  element 
approach.  Only  linear  elastic  properties  were  examined. 
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Aboudi[2lJ  developed  a  theoretical  model  for  an  imperfect  bond  based 
on  the  simplifying  assumptions  of  Jones  and  Whittier[22],  who  expressed 
an  imperfect  bond  as  a  thin  elastic  film  containing  two  parameters.  Each 
parameter  represents  the  stiffness  in  the  tangential  and  normal  directions, 
respectively.  Aboudi’s  results  agreed  favorably  with  those  of  Benveniste 
[23;  24]  who  used  a  concentric  cylinder  three  phase  model.  Other  theoretical 
models  of  debonding  have  appeared  for  determining  a  weak  bond’s  affect 
on  longitudinal  shear{25],  for  analyzing  wave  propagation  in  unidirectional 
composites  [26],  and  for  particulate  composites  possessing  imperfectly 
bonded  spheres[27]  have  also  appeared. 

In  addition,  Pagano  and  Tandon[28-31]  used  a  circular  concentric  cyl¬ 
inder  elastic  model  they  had  previously  developed  similar  to  the  self  con¬ 
sistent  model  to  analyze  interfacial  debonding.  Complete  separation  was 
modeled  by  replacing  the  fibers  with  cylindrical  voids,  and  a  slipping 
condition  was  simulated  by  forcing  the  shear  traction  at  the  interface  to 
vanish. 

All  of  these  linear  elastic  methods  meet  with  the  same  difficulty. 

When  attempting  to  model  an  imperfect  bond  using  linear  elastic  methods 
alone,  a  weak  interfacial  bond  condition  necessarily  exists  in  compression 
as  well  as  in  tension  since  the  stiffness  must  remain  constant.  In  other 
words,  such  an  analysis  will  result  in  the  interface  being  incapable  of  sus¬ 
taining  a  compressive  load,  and  consequently,  the  fiber  and  matrix  will  be 
allowed  to  pass  through  one  another.  Therefore,  no  constraints  will  exist 


during  the  application  of  thermal  attains,  and  hence,  no  thermal  stresses 


will  result.  This  is  a  major  drawback  to  the  models  discussed  so  far  when 
analyzing  high  temperature  metal  matrix  composites. 
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In  an  effort  to  model  thermal  residual  stresses  Hopkins  and  Chamis[32] 
developed  a  set  of  micromechanics  eolations  for  high  temperature  metal 
matrix  composites  based  on  a  concentric  cylinder  mechanics  of  materials 
approach.  However,  in  order  to  include  thermal  effects,  it  was  necessary 
to  assume  a  strong  bond  between  the  constituents. 

Including  both  a  weak  fiber/matrix  bond  and  thermal  residual  stresses 
in  a  micromechanical  analysis  is  necessarily  a  nonlinear  problem.  Resid¬ 
ual  compressive  stresses  at  the  interface  must  be  overcome  before  the  ma¬ 
trix  will  separate  from  the  fiber.  As  mentioned  in  the  first  chapter,  this 
was  experimentally  observed  by  Johnson  et.  al.[l],  A  distinct  “knee”  in 
the  transverse  load-displacement  curve  occurred.  It  was  postulated  that  a 
weaker  modulus  existed  after  fiber/matrix  separation  resulting  in  the  ap¬ 
proximate  bilinear  behavior. 

Residual  thermal  stresses  are  not  the  only  mechanism  controlling 
fiber/matrix  separation.  The  fiber/matrix  bond  seems  to  possess  a  finite 
ultimate  strength.  Naik  et.  al .  [33]  found  that  during  initial  loading  of  an 
actual  specimen  the  characteristic  “knee”  occurred  at  a  slightly  higher 
load  than  in  subsequent  cycles.  The  additional  loading  before 
nonlinearity  occurred  during  the  first  cycle  was  assumed  to  be  the  load  re¬ 
quired  to  fail  the  fiber/matrix  bond.  Since  bond  failure  had  already  oc¬ 
curred,  subsequent  loadings  demonstrated  bilinear  behavior  earlier.  Other 
than  for  purposes  of  discussion,  a  finite  ultimate  strength  of  the  fiber/ma¬ 
trix  bond  will  not  be  considered  in  the  present  study.  Since  it  is  present 
during  initial  loading  only,  it  will  not  be  of  concern  in  most  structural  ap¬ 
plications. 
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Highsmith  and  Naik[34]  attempted  to  determine  the  onset  of  separation 
analytically  using  both  a  concentric  cylinder  self  consistent  model  and  a 
linear  elastic  finite  element  analysis.  Again,  due  to  the  linear  elastic  na¬ 
ture  of  the  analysis  it  was  necessary  to  assume  a  strong  bond  for  the  cal¬ 
culations,  and  as  a  result  the  analysis  could  not  determine  the  characteris¬ 
tics  beyond  the  point  of  separation.  If  one  desires  to  determine  the  com¬ 
posite  characteristics  after  this  point,  as  a  minimum,  a  bilinear  capability 
must  be  present  :.t  the  interface.  Therefore,  using  the  types  of  linear  anal¬ 
yses  discussed  so  far  is  not  practical. 

However,  even  a  bilinear  elastic  analysis  will  not  allow  for  energy 
dissipation,  and  it  is  apparent  that  there  is  a  significant  amount  of  energy 
dissipation  during  transverse  loading  and  unloading  of  a  titanium-based 
composite.  If  only  a  small  portion  of  the  matrix  yields  during  loading, 
the  unloading  phenomena  will  be  dominated  by  the  elastic  stresses  in  the 
matrix  which  contains  significant  thermal  residual  strain.  Therefore,  upon 
unloading,  the  elastic  stresses  in  the  matrix  will  tend  to  close  the  separat¬ 
ed  fiber/matrix  interface  before  complete  unloading  occurs.  However,  a 
significant  amount  of  energy  will  be  lost  at  the  interface  due  to  frictional 
effects,  crack  propagation,  and  plastic  flow[l]. 

Nimmer[2]  performed  a  nonlinear  analysis  of  a  weak  interfacial  bond 
coupled  with  thermal  residual  stresses  using  both  a  simplified  square  fiber 
model  and  a  finite  element  analysis.  He  examined  a  monotonically  in¬ 
creasing  transverse  normal  load  and  compared  the  characteristics  of  a 
Stronp  and  weak  interfarial  bond.  A  uniform  sanare  arrav  of  fibers  was 
assumed  and  then  symmetry  used  to  isolate  one  quarter  of  a  unit  fiber/ma¬ 
trix  cell.  The  fiber  was  assumed  to  be  a  circular  cylinder  with  no  voids 
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between  it  and  the  surrounding  matrix.  His  work  agreed  at  least  qualita¬ 
tively  with  the  experimental  results  of  Johnson  et.  al .  [1] .  A  key  part  of  the 
analysis  was  the  incorporation  of  a  contact  element  at  the  fiber/matrix  in¬ 
terface.  This  type  of  two  noded  element  allows  for  high  stiffness  in  com¬ 
pression  and  extremely  low  stiffness  in  tension.  Also,  a  friction  coeffi¬ 
cient  may  be  included  for  transverse  slippage  along  the  interface. 

Figure  1  presents  Nimmer’s  results  on  the  effect  of  thermal  residual 
stresses  in  the  presence  of  a  weak  fiber/matrix  bond.  A  35%  fiber  volume 
fraction  was  used,  and  the  zero  stress  state  was  taken  to  be  1000°C.  The 
model  was  subsequently  cooled  to  23°C  followed  by  a  monotonically  in¬ 
creasing  load.  For  the  weak  interface  solution  the  Coulomb  friction 
coefficient(friction  coefficient  between  two  dry  surfaces)  was  taken  to  be 
1.0,  and  the  material  properties  are  given  in  Table  1.  The  strong  interface 
model  does  not  experience  yielding  until  approximately  950  MPa  while  the 
weak  interface  solution  shows  an  apparent  yielding  at  less  than  200 MPa. 
Nimmer  demonstrated  that  this  point  of  apparent  yield  was  actually 
fiber/matrix  separation.  True  yielding  of  the  matrix  in  the  weak  interface 
model  did  not  occur  until  well  above  200 MPa.  Thus,  the  weak  model  was 
still  elastic  well  after  it  experienced  the  nonlinearity. 

Subsequently,  Nimmer[3]  attempted  a  more  in-depth  investigation  of 
his  finite  element  model  and  compared  it  to  experimental  data(Figure  2). 

A  composite  of  silicon  carbide  fiber  and  titanium  aluminide  matrix  was 
used  for  the  experiment,  and  its  associated  constituent  properties  were 
input  into  the  analysis.  An  attempt  was  made  in  the  computations  to 
match  the  experimental  specimen’s  fiber  distribution  aspect  ratio  and  vol¬ 
ume  fraction.  Finally,  analysis  and  experiment  were  brought  through  a 
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type  of  interface: 

- strong 

- weak 


0  0.001  0.002  0.003  0.004  0.005  0.006  0.007  0.008 

Strain 

Figure  1.  Nimmer’s  Results  for  Strong  vs.  Weak  Interface,  Plane  Stress 
Finite  Element  Solution (2] 


Table  1.  Material  Properties  Used  in  the  Finite  Element  Solutions  of 
Figure  1. 


E  (GPa) 

V 

a  (“C)-1 

Oys  (MPa) 

Fiber 

400. 

0.3 

5.0x10-6 

- 

Matrix 

100. 

0.3 

lO.OxlO-6 

1000. 

Stress 

(MPa) 


Figure  2.  Nimmer’s  Results  of  Analysis  vs.  Experiment,  3-D  General¬ 
ized  Plane  Strain  Finite  Element  Solution[3] 


similar  thermomechanical  loading  sequence  consisting  of  cooldown  fol¬ 
lowed  by  loading,  unloading,  and  then  reloading.  The  experimental  re¬ 
sults  exhibited  almost  no  permanent  strain  upon  unloading  confirming  that 
little  matrix  plasticity  occurs.  In  addition,  the  small  finite  strength  of  the 
fiber/matrix  bond  is  demonstrated  in  the  difference  between  the  initial 
loading  curve  and  reload. 

Nimmer’s  generalized  plane  strain  model  does  an  excellent  job  of 


predicting  the  initial  modulus  as  well  as  the  point  of  fiber/matrix  separa¬ 


tion.  Since  the  experimental  specimen  demonstrates  interfacial  bond  fail¬ 
ure  upon  initial  loading  and  the  analysis  models  the  interface  as  weak,  the 
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initial  loading  in  the  analysis  should  compare  most  favorably  with  the  re¬ 
loading  of  experiment.  Nevertheless,  the  secondary  stiffness  after  separa¬ 
tion  is  considerably  weaker  than  the  empirical  data,  and  in  addition,  upon 
unloading  the  model  fails  to  return  to  its  approximate  original  strain  state 
as  observed  experimentally.  This  indicates  that  modeling  the  fiber/matrix 
contact  as  perfectly  weak  is  not  an  adequate  representation  of  the  actual 
mechanism.  Hence,  there  must  be  additional  constraints  existing  at  the 
interface  after  separation.  Also,  since  the  model  exhibits  significant  plas¬ 
tic  strain  after  unloading,  it  does  not  correctly  model  a  matrix  dominated 
by  elastic  forces. 

The  only  means  of  energy  dissipation  in  Nimmer’s  model  were  trans¬ 
verse  friction  along  the  fiber/matrix  contact  and  plastic  flow  of  the  ma¬ 
trix.  However,  the  experimental  data  indicates  that  a  significant  amount 
of  energy  dissipation  occurs  even  when  the  matrix  is  dominated  by  elastic 
forces.  Since  very  little  of  the  matrix  yields,  the  bulk  of  the  energy  loss 
must  occur  at  the  fiber/matrix  interface.  It  is  essential  to  model  this  ener¬ 
gy  loss  if  correct  load  histories  of  the  composite  are  to  be  reproduced  in  a 
micromechanical  analysis. 

Energy  absorption  is  a  prime  factor  in  determining  the  life  cycle  of  a 
material.  In  many  cases,  fatigue  failure  criteria  is  based  on  the  energy 
dissipated  during  loading  and  unloading [35].  Therefore,  before  reliable 
life  prediction  techniques  for  high  temperature  metal  matrix  composites 
may  be  obtained,  a  strong  understanding  of  the  energy  dissipation  process 
is  needed.  To  develop  accurate  models,  parametric  studies  must  first  be 
performed  using  various  interface  conditions  to  determine  what  mecha¬ 
nisms  must  be  included. 
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The  results  of  such  studies  will  demonstrate  the  effect  of  various  in¬ 
terface  conditions  on  the  composite  characteristics.  In  the  study  present¬ 
ed  here  an  attempt  will  be  made  to  determine  if  the  energy  dissipation 
mechanisms  existent  at  the  interface  may  be  adequately  modeled  by  repre¬ 
senting  the  reaction  zone  as  a  material  with  different  material  properties 
as  compared  to  the  matrix. 

Much  work  remains  to  be  done  in  understanding  high  temperature 
metal  matrix  composites.  In  this  effort,  micromechanical  analysis  can 
provide  both  the  ability  to  predict  macroscopic  characteristics  of  the  com¬ 
posite  and  furnish  a  clearer  understanding  into  the  mechanisms  surround¬ 
ing  interfacial  debonding.  Once  equipped  with  a  stronger  understanding 
of  these  microscopic  mechanisms,  micromechanical  analysis  can  bridge 
the  gap  between  experimental  data  and  developing  reliable  life  cycle  pre¬ 
diction  capabilities.  Therefore,  it  is  imperative  that  a  reasonable  model  of 
the  interface  zone  be  developed. 
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Ill  -  Methods  of  Analysis 


3.1  Hardware  and  Software 

Two  types  of  analysis  will  be  used  in  this  study:  the  finite  element 
method  and  a  simplified  one  dimensional  model.  A  nonlinear  material  fi¬ 
nite  element  method  will  receive  the  most  attention.  However,  the  simpli¬ 
fied  one  dimensional  model  will  also  be  employed  for  rapid  calculation  of 
trends  in  the  data.  Since  this  will  be  a  parametric  study  where  many  con¬ 
ditions  will  be  examined,  using  a  finite  element  method  alone  would  be 
quite  tedious.  Each  nonlinear  solution  may  take  anywhere  from  2  to  over 
24  hours  of  CPU  time.  When  one  also  considers  the  time  required  for  grid 
generation  and  application  of  constraints,  it  becomes  a  task  far  from  at¬ 
tractive.  Additionally,  unless  another  means  of  verification  is  available, 
using  a  simplified  analysis  alone  makes  one  to  question  the  results. 
Therefore,  a  mixture  of  finite  element  and  one  dimensional  solutions  will 
be  used. 

Computer  system  hardware  included  the  aeronautical  engineering 
department’s  SUN  4/330  workstation  for  the  finite  element  calculations 
and  pre/post  processing  of  results.  When  the  workstation  was  unavail¬ 
able,  the  computer  directorate’s  VAX  11/785  (Hercules)  was  employed. 
Finally,  an  Apple  Macintosh  IIx  was  used  for  the  one  dimensional  model 
calculations. 

Important  software  consisted  of  the  Macneal-Schwendler  Corpora¬ 
tion’s  MSC/NASTRAN  finite  element  analysis  program  and  the  Structural 
Dynamics  Research  Corporation’s  SDRC-IDEAS  computer  aided  de- 
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sign/engineering  package.  The  SDRC-IDEAS  package  was  used  only  as  a 
pre/post  processor  for  MSC/NASTRAN.  Additionally,  the  one  dimension¬ 
al  calculations  were  performed  using  Mathematic  a™,  a  versatile  mathe¬ 
matics  program  capable  of  symbolic  operations. 

The  remainder  of  this  section  will  present  a  brief  discussion  of  the  ca¬ 
pabilities  of  MSC/NASTRAN  for  nonlinear  material  analysis,  software 
utilization,  and  some  of  the  difficulties  encountered.  MSC/NASTRAN’s 
nonlinear  material  static  analysis  solver(solution  sequence  66)  uses  a 
modified  Newton-Raphson  approach  for  convergence.  The  tangent  method 
is  employed,  and  the  user  maintains  the  option  of  when  to  update  the  stiff¬ 
ness  matrices.  Other  convergence  parameters  available  to  the  user  include 
the  number  of  load  increments  per  subcase,  the  type  of  convergence 
test(displacement,  load,  or  energy  error  fractions),  and  the  error  tolerance 

limits.  Both  load  error  fraction,  ep,  and  energy  error  fraction,  ew,  were 
incorporated  to  control  convergence.  These  quantities  are  defined  in  the 
MSC/NASTRAN  Application  Manual  [36:  2.14-52].  The  solution  was  as¬ 
sumed  to  have  converged  when  ep^1.0xl0'3  and  ew^1.0x!0-7. 

In  linear  elastic  analysis  the  loading  sequence  is  immaterial,  but  for 
nonlinear  elastic-plastic  problems  the  solution  is  dependent  on  the  load 
history.  NASTRAN  incorporates  a  load  history  by  applying  the  load 
through  user  specified  subcases.  Each  subcase  may  further  be  increment¬ 
ed  for  more  stable  convergence.  Once  a  solution  for  the  subcase  is 
achieved  the  next  subcase’s  calculations  are  begun  starting  from  the  pre¬ 
vious  solution. 

Several  types  of  yield  criteria  are  available  in  MSC/NASTRAN.  The 
Von  Mises  criterion  gives  good  results  for  most  metals,  and  since  a  metal 
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matrix  material  was  to  be  analyzed,  this  criterion  was  chosen  for  the 
present  study.  The  Von  Mises  stress  is  defined  as  follows: 

0VON  MISES  =  /\J1^  (Ox  -  CJyj^  +  (CJy  -  Gz)2  +  (Gz  -  Ojf"  +  6  (ixy2  +  Ty z2  +  T^2)] 

(1) 

where  Gx>  Gy,  and  Gz  are  the  normal  stresses  along  any  orthogonal  set  of 
axes,  and  Txy,  Tyz,  and  Tzx  are  the  respective  shear  stresses  for  the  axes. 
The  Von  Mises  criterion  states  that  yielding  occurs  when  the  Von  Mises 
stress  exceeds  the  yield  stress. 

MSC/NASTRAN  incorporates  the  Prandtl-Reuss  flow  relations,  and 
possesses  the  capability  for  either  isotropic,  kinematic,  or  a  combined 
hardening.  Since  a  loading  and  unloading  sequence  is  planned,  kinematic 
hardening  was  chosen  in  order  to  account  for  the  Bauschinger  effect  in  the 
simplest  manner. 

The  types  of  elements  used  in  this  study  included  a  three  noded  con¬ 
stant  strain,  constant  curvature  triangular  shell  element,  a  four  noded  iso¬ 
parametric  shell  element,  a  six  noded  three  dimensional  isoparametric 
wedge  element,  an  eight  noded  three  dimensional  isoparametric  brick  ele¬ 
ment,  and  a  two  noded  gap  element.  The  user  input  for  the  shell  and  solid 
element’s  material  properties  is  the  one-dimensional  stress-strain  curve 
for  tension.  The  gap  element  provides  a  unique  capability  necessary  for 
the  present  study.  It  is  an  element  intended  to  model  surfaces  which  may 
come  into  contact.  It  connects  two  grid  points(nodes)  which  may  be  ini¬ 
tially  coincident,  and  the  local  element  coordinate  system  defines  a  con¬ 
tact  plane  which  does  not  change  with  deflection.  When  compressive 
normal  stress  exists,  the  gap  element  may  possess  a  longitudinal  stiffness 
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comparable  to  the  surrounding  material  and  carry  a  transverse  shear  load 
of  a  magnitude  less  than  the  normal  load  times  the  coefficient  of  friction. 
When  the  element  is  in  tension,  its  longitudinal  stiffness  may  be  dramati¬ 
cally  reduced  to  approximately  zero  and  it  is  no  longer  capable  of 
resisting  a  transverse  shear  load. 

A  good  portion  of  the  grid  generation  and  data  reduction  was  pei- 
formed  using  the  SDRC-IDEAS  software.  Using  this  program  it  was  pos¬ 
sible  to  minimize  the  number  oi  manual  inputs  to  a  NASTRAN  data  deck. 
A  few  minor  difficulties  arose.  It  was  found  that  IDEAS  was  unable  to 
write  out  temperature  loads,  and  incorrect  Direct  Matrix  Abstraction  Pro¬ 
gramming  (DMAP)  statements  for  creating  the  necessary  post  processing 
files  from  NASTRAN  were  listed  in  the  IDEAS  manual  [38].  For  the  cor¬ 
rect  DMAP  statements  see  appendix-A  which  lists  a  sample  MSC/NAS- 
TRAN  input  deck  used  in  this  study. 

A  couple  of  significant  software  problems  in  using  MSC/NASTRAN 
are  worth  mentioning  for  future  reference.  Rigid  body  elements  are  ig¬ 
nored  when  thermal  loads  are  present.  Therefore,  rather  than  using  rigid 
body  elements  to  ensure  constant  displacement,  multipoint  constraints 
(MPCs)  must  be  written.  Also,  in  the  nonlinear  material  solution  se¬ 
quence  there  was  a  major  discrepancy  in  the  solutions  between  MSC/NAS¬ 
TRAN  version  65C  and  the  recently  released  ve  *'4on  66A.  Since  the  65C 
results  demonstrated  excellent  agreement  with  Nimmer’s  solutions  calcu¬ 
lated  using  the  finite  element  program  ADJ.NA,  all  recults  presented  in 
this  study  are  from  the  65C  version.  Due  to  this  discovery  the  MacNeal- 
Schwendler  Corporation  is  presently  looking  into  this  discrepancy,  and 
until  the  problem  is  overcome  the  author  recommends  using  version  65C 
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for  any  nonlinear  material  solutions. 

3.2  General  Information  on  Models  Used  in  this  Study 

In  order  to  perform  the  investigation  several  assumptions  were  neces¬ 
sary.  First,  a  regular  rectangular  array  of  fibers  was  assumed  for  all  mod¬ 
els.  A  unit  fiber/matrix  cell  or  representative  volume  element  (RVE)  con¬ 
sists  of  a  single  fiber  and  the  surrounding  matrix  material.  Symmetry  was 
used  in  the  analysis  to  isolate  one  quarter  of  a  unit  fiber/matrix  cell  for 
the  calculations  (Figure  3).  The  fiber  diameter  was  taken  to  be  0.14  mm, 
and  the  remaining  dimensions  of  ihe  analysis  cell  were  determined  from 
the  fiber  volume  fraction  and  aspect  ratio.  The  material  properties  of  the 
fiber,  matrix,  and  interphase  zone  were  assumed  to  be  isotropic. 

Boundary  conditions  along  the  external  faces  of  the  analysis  cell  were 
chosen  as  the  simplest  to  satisfy  compatibility.  Thus,  -he  lower  face  is 
fixed  in  the  vertical  direction  only,  and  the  left  face  in  the  horizontal  di¬ 
rection  only.  The  other  two  in-plane  faces  are  free  to  move  parallel  to 
their  faces,  but  perpendicularly  they  are  constrained  to  move  uniformly. 
The  loading  for  all  models  consisted  of  a  cooldown  from  a  specified  pro¬ 
cessing  temperature  to  room  temperature  followed  by  a  transverse  normal 
load  applied  to  the  right-hand  face. 

Two  representative  volume  elements  will  be  used  throughout  this 
study:  a  35%  fiber  volume  fraction  with  an  assumed  square  fiber 
array(aspect  ratio=1.0),  and  a  32%  fiber  volume  fraction  of  a  rectangular 
array  with  an  aspect  ratio  of  1.2.  The  35%  fiber  volume  fraction  square 
array  representative  volume  element  which  at  times  will  be  referred  to  as 
RVE#1  was  used  as  a  generic  model  for  determining  how  the  composite 


Figure  3.  Basic  Cell  Used  in  Analysis 


behaved  under  various  interface  conditions.  The  properties  used  were 
similar  to  those  used  previously  by  Nimmer[2].  Since  this  is  a  parametric 
study,  RV£#1  will  be  employed  most  often.  The  32%  fiber  volume  frac¬ 
tion  rectangular  array  representative  volume  element  with  the  1.2  aspect 
ratio(RVE#2)  is  for  specific  comparison  with  an  actual  specimen  from 
previously  published  experimental  data  where  micrographs  were  used  to 


obtain  the  correct  fiber  volume  fraction  and  aspect  ratio[3]. 

3.3  Finite  Element  Models  for  Strong  and  Weak  Interfaces 

The  strong  and  weak  interface  models  were  created  mainly  for  com¬ 
parison  with  Nimmer’s  previous  work  as  a  check  for  the  calculations  and 
to  provide  a  baseline.  However,  in  a  few  of  the  cases  there  was  some  ex¬ 
trapolation  on  his  work,  and  in  one  instance  an  additional  constraint  ap¬ 
plied.  The  strong  and  weak  interface  finite  element  models  consisted  of 
both  two  dimensional  and  three  dimensional  grids.  Both  plane  stress  and 
plane  strain  solutions  were  performed  on  the  two  dimensional  models. 

For  the  strong  bond  calculations  only  RVE#1  was  used.  The  two  di¬ 
mensional  strong  interface  model  consisted  of  quadrilateral  isoparametric 
four  noded  shell  elements  and  triangular  constant  strain  constant  curva¬ 
ture  three  noded  shell  elements.  Since  only  membrane  effects  were  need¬ 
ed,  all  rotations  and  out  of  plane  displacements  were  removed  through 
single  point  constraints.  Also,  since  the  strong  interface  model  experienc¬ 
es  little  or  no  nonlinearities  at  the  load  levels  considered(no  fiber/matrix 
separation  occurs  and  the  transverse  load  is  less  than  half  of  the  matrix 
yield  stress),  then  only  a  linear  analysis  was  necessary.  The  finite  ele¬ 
ment  grid  is  shown  in  Figure  4  and  was  biased  towards  regions  of  expect¬ 
ed  maximum  stress.  This  biasing  ratio  is  approximately  3  to  1. 

The  boundary  conditions  along  the  four  faces  were  applied  through 
single  point  and  multipoint  constraints  as  well  as  nodal  applied  loads. 
Single  point  constraints  are  simply  specified  displacements  at  particular 
nodes,  while  a  multipoint  constraint  is  a  linear  equation  that  relates  a  set 
of  degrees  of  freedom  to  each  other.  Thus,  a  multipoint  constraint  for  the 
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Figure  4.  35%  Fiber  Volume  Fraction  Square  Array  Finite  Element 
Grids  for  Both  Strong  and  Weak  Interfacial  Bonds 
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set  of  degrees  of  freedom,  {u1,u2,...,un} ,  will  take  the  form 


n 

X  aiui  =  0 

i  =  1 


(2) 


where  aj  are  constant  coefficients. 

The  nodes  on  the  lower  face  were  restrained  in  the  vertical  direction 
only,  and  the  nodes  on  the  left  face  were  restrained  in  the  horizontal  di¬ 
rection  only.  Multipoint  constraints  were  applied  on  the  right  and  upper 
fa-.es  to  provide  for  equal  displacement.  All  nodes  on  the  right  face  are 
tied  through  a  multipoint  constraint  to  the  node  on  the  lower  right  corner 
in  such  a  way  that  they  will  displace  equally  in  the  horizontal  direction. 
The  same  is  true  for  the  upper  face  in  the  vertical  direction  where  the 
node  on  the  upper  left  corner  is  used.  Transverse  loading  was  applied  at 
the  lower  right-hand  node.  Determining  equivalent  nodal  loads  for  an 
evenly  applied  face  pressure  was  not  necessary  since  all  the  nodes  along 
the  face  will  displace  equally. 

To  create  the  finite  element  grio  for  the  three-dimensional  strong  bond 
model  the  two  dimensional  grid  was  simply  extruded  into  three  dimen¬ 
sional  space  by  a  thickness  of  .01  mm.  Therefore,  it  was  made  up  of  eight 
noded  brick  elements  and  six  noded  wedge  elements.  The  purpose  of  the 
three  dimensional  model  was  to  analyze  behavior  under  a  generalized 
plane  strain  condition.  Since  the  composite  longitudinal  direction  is  the 
out-of-plane  direction  in  the  analytical  models,  it  is  assumed  that  this  di¬ 
rection  is  very  thick,  and  hence,  a  generalized  plane  strain  solution  is  the 
most  physically  realistic.  No  change  of  the  in-plane  boundary  conditions 
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were  necessary  for  the  three-dimensional  model,  but  an  additional  degree 
of  freedom  in  the  out-of-plane,  z-direction,  must  be  dealt  with.  These  ad¬ 
ditional  boundary  conditions  consist  of  the  back  face  being  fixed  in  the  z- 
direction  and  the  nodes  on  the  front  face  being  constrained  to  displace 
uniformly  in  the  z-direction  through  the  use  of  multipoint  constraints.  All 
the  strong  bond  calculations  provided  a  good  reference  for  the  weak  inter¬ 
face  solutions. 

The  weak  interface  models  contained  finite  element  grids  of  both 
RVE#1  and  RVE#2.  The  finite  element  grid  for  the  weak  bond  RVE#1 
model  was  identical  to  the  strong  bond  model  except  for  the  incorporation 
of  gap  elements  to  allow  for  the  weak  fiber/matrix  interface  condition. 

The  weak  model  possessed  coincident  nodes  at  the  interface  where  one  set 
is  attached  to  the  matrix  and  the  other  to  the  fiber.  These  coincident  in¬ 
terface  nodes  are  then  connected  to  one  another  through  gap  elements  that 
prevent  the  fiber  and  matrix  from  passing  through  one  another  but  on  thf 
other  hand  offer  no  resistance  under  separation.  In  addition,  a  wei.*  -r 
face  model  of  RVE#2  was  also  developed.  Its  finite  element  gr’d  ■ 
played  in  Figure  5,  and  except  for  the  dimensions  all  other  characteristics 
of  it  are  similar  to  the  RVE#1  model.  Finally,  as  would  be  expected,  the 
boundary  conditions  for  the  weak  bond  models  were  identical  to  those  dis¬ 
cussed  in  the  previous  paragraphs  for  the  strong  interface  models. 

In  order  to  distinguish  the  effect  of  interfacial  slippage  between  the 
fiber  and  matrix,  an  additional  constraint  which  prevents  any  transverse 
slippage  along  the  interface  but  still  allows  for  separation  was  developed 
for  the  RVE#2  model.  Again  multipoint  constraints  were  utilized  to  only 
allow  the  coincident  nodes  to  displace  perpendicularly  from  the  interface. 
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Figure  5.  32%  Fiber  Volume  Fraction  Rectangular  Array(Aspect 

Ratio=1.2)  Finite  Element  Grids  for  Weak  Interfacial  Bond 
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This  constraint  could  be  optionally  applied,  so  it  was  possible  to  make  di¬ 
rect  comparisons  between  the  free  slip  and  no  slip  conditions. 

3.4  Models  for  Incorporating  an  Interphase  Zone 

A  considerable  portion  of  this  study  was  an  attempt  to  model  the  im¬ 
perfections  that  exist  in  the  reaction  zone  by  incorporating  a  third  zone  at 
the  fiber/matrix  interface  possessing  isotropic  material  properties.  This 
zone  will  be  modeled  with  various  material  properties  in  order  to  simulate 
the  imperfect  bond  characteristics  and  to  allow  for  energy  dissipation  in  a 
loading  and  unloading  cycle.  Both  a  simplified  one  dimensional  model 
and  finite  element  analysis  will  be  used  for  analyzing  the  effect  of  an  in¬ 
terphase  zone. 

The  one  dimensional  model  is  patterned  after  a  bilinear  elastic  analy¬ 
sis  proposed  by  Nimmer  for  a  weak  bond[2].  However,  a  higher  level  of 
sophistication  is  needed  to  model  an  interphase  zone  and  its  elastic-plastic 
characteristics.  A  visualization  of  the  model  is  depicted  in  Figure  6 
where  a  square  fiber  is  assumed  and  the  entire  model  is  divided  into  five 
regions  consisting  of  two  matrix  regions,  two  interphase  zone  regions,  and 
a  single  fiber  region.  The  stress  and  strain  is  assumed  constant  through¬ 
out  each  region,  and  to  reduce  the  equations  into  one  dimension  both  the 
shear  and  Poisson’s  effect  are  neglected.  Hence,  only  the  stress  and  strain 
in  the  direction  of  applied  load  is  necessary  which  in  this  case  is  the  x-di- 
rection.  Since  the  left  face  is  fixed,  the  strains  in  each  region  may  be  re¬ 
lated  to  the  displacements  shown  in  Figure  6  as  follows: 
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(4) 


eii  = 


Ull-UF 

b 


_  „  UM1-UI1 

Mj  c 


(5) 


(6) 


EM2  " 


UM2 

(a  +  b  f  c) 


(7) 


where  u  is  displacement,  E  is  strain,  and  subscripts  are  the  specified  re¬ 
gions  shown  in  Figure  6. 

Both  the  fiber  and  matrix  regions  are  assumed  to  be  linearly  elastic, 
and  therefore  the  only  material  properties  needed  for  these  regions  are 
Young’s  modulus  and  thermal  coefficient  of  expansion.  The  interphase 
zone  is  assumed  to  be  elastic-plastic  with  isotropic  hardening  during  plas¬ 
tic  deformation.  Therefore,  the  material  properties  required  for  the  inter¬ 
phase  zone  are  Young’s  modulus,  thermal  coefficient  of  expansion,  initial 
yield  stress,  and  a  strain  hardening  parameter.  Incorporating  the  plastic 
strains  in  both  interphase  zone  regions  and  denoting  the  cooldown  temper¬ 
ature  by  T  and  the  processing  temperature  by  Tp  the  one  dimensional 
stress-strain  relations  for  the  model  become: 

eF=|E+aF(T.Tp)  (g) 

Ell  =  EEi!  +  EPll  (9) 
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(10) 


EEl,-^  +  Ctl(T-Tp) 

£Ml=^L+ClM(T'Tp) 

£i2  =  £ei2  +  £pi2 

EEl2  =  ^+  CCI(T-Tp) 
£M2  =  +  ttM  (T  ' Tp) 


(11) 

(12) 


(13) 


(14) 


where  G  is  stress,  a  is  the  thermal  coefficient  of  expansion,  E  is  Young’s 
modulus,  the  subscripts  E  and  P  denote  elastic  and  plastic  strains,  and  the 
other  symbols  and  subscripts  are  the  same  as  in  Eqs  (3-7). 

Designating  the  applied  stress  on  the  right-hand  face  by  Oa,  the  equi¬ 
librium  equations  for  this  simplified  model  are: 


Of  a  +  Oi2  b  =  Oij  (a  +  b) 

(15) 

Oij  (a  +  b)  =  Omi  (a  +  b) 

(16) 

OMi  (a  +  b)  +  0  m  2  c  =  Oa(a  +  b  +  c) 

(17) 

There  are  nineteen  unknowns  in  this  set  of  fifteen  equations,  so  addi¬ 
tional  information  is  required.  Two  more  equations  may  be  obtained  by 
applying  certain  constraints  on  displacements.  These  constraints  come 


33 


from  the  boundary  condition  on  the  right-hand  face  and  compatibility  of 
the  interior.  They  are: 


Ui2  =  Up 

oo 

uM2  =  uMi 

(19) 

The  problem  is  still  in  need  of  more  information.  This  comes  in  the 
form  of  the  one  dimensional  yield  criteria  and  the  strain  hardening  rela¬ 
tion.  The  yield  criteria  states  that  if  the  interface  stress  exceeds  the  yield 
stress,  then  plastic  deformation  occurs,  otherwise  the  material  behaves 
elastically.  Since  plastic  solutions  are  dependent  on  the  load  history,  the 
load  must  be  applied  in  incremental  steps  so  that  the  plastic  strains  from 
the  previous  load  step  may  be  used  in  the  present  calculation.  The  one  di¬ 
mensional  model  presented  here  employs  the  following  solution  sequence. 
If  yielding  does  not  occur  during  an  applied  load,  then  the  previous  load 
step’s  plastic  strains  are  used  in  the  above  17  equations  to  solve  for  the 
17  unknowns.  If  yielding  does  occur,  then  the  yielded  region’s  stress  is 
set  equal  to  the  yield  stress  of  the  region  which  could  be  either  in  com¬ 
pression  or  tension,  and  its  associated  plastic  strain  becomes  the  unknown 
to  be  solved  for.  The  yield  stress  for  each  interphase  region  is  given  by 
the  following  strain  hardening  relations. 

aysij=  aysi+ H' £Pll  (20) 

Oysi2=  oysi+ H' epl2  (2i) 

where  Cfysj  is  the  initial  yield  stress  for  the  interphase  material,  and  H'  is 
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the  strain  hardening  parameter. 

A  Mathematica™  program  was  written  to  solve  Eqs  (3-21)  for  all 
possible  interphase  zone  yielding  combinations  of  the  two  regions  both  in 
compression  and  tension.  Any  transverse  rmal  loading  sequence  may  be 
applied  to  the  mathematical  model,  but  in  this  study  only  a  single  loading 
and  unloading  sequence  will  be  used.  The  program  is  listed  in  Appendix- 
13.  During  an  applied  load  step  the  stresses  are  initially  assumed  to  be 
elastic,  and  upon  evaluation  of  the  stresses  in  the  interphase  regions  a 
judgment  is  made  through  the  use  of  a  logic  sequence  as  to  what  set  of 
equations  apply.  This  procedure  is  continued  until  two  successive  solu¬ 
tions  are  identical. 

The  one  dimensional  model  provided  a  good  supplement  to  the  finite 
element  solutions  by  providing  rapid  calculations  of  how  varying  a  specif¬ 
ic  property  or  parameter  generally  affected  the  overall  structure.  As  a  re¬ 
sult,  more  effective  use  of  the  finite  element  solutions  were  possible.  In 
addition,  it  provided  a  good  check  for  the  finite  element  calculations  by 
preventing  gross  errors  from  being  overlooked  as  might  occur  if  no  other 
analysis  tool  was  available. 

Similar  to  the  perfectly  weak  bond  calculations,  the  finite  element  so¬ 
lutions  that  incorporated  an  interphase  zone  contained  both  two  and  three 
dimensional  models.  The  same  boundary  conditions  were  employed,  and 
the  three  dimensional  models  were  again  obtained  by  extruding  the  two 
dimensional  ones  into  three  dimensional  space.  Since  the  effect  of  the  in¬ 
terphase  zone  size  was  to  be  examined,  separate  grids  had  to  be  created 
for  each  size  investigated.  Three  sizes  were  examined.  Their  thicknesses 
were  5%,  7.5%,  and  10%  of  the  fiber  radius.  Such  thicknesses  were  cho- 
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Figure  7.  Coarse  and  Fine  Grids  for  Interphase  Zone  Size  5%  of  Fiber 
Radius 
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Figure  8.  Finite  Element  Grids  for  Interphase  Zone  Sizes  of  7.5%  and 
10%  of  Fiber  Radius 
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3.5 


Combined  lnterohase  Zone  and  Weak  Bond  Models 


All  the  models  to  this  point  have  assumed  that  the  fiber/matrix  inter¬ 
facial  region  possesses  identical  properties  at  all  points  around  the  fiber. 
Most  likely  this  is  not  the  case  in  an  actual  specimen.  Portions  of  the  in¬ 
terface  may  represent  a  perfectly  weak  bond  while  other  portions  may  still 
offer  resistance  either  in  separation  or  an  irregular  fracture  surface  may 
arrest  transverse  slippage  at  the  interface.  Therefore,  additional  models 
were  developed  to  analyze  a  combination  of  fiber/matrix  interface  condi¬ 
tions.  These  models  consist  of  a  section  of  the  interface  taken  to  be  per¬ 
fectly  weak  while  the  remainder  is  assumed  to  behave  as  a  third  phase  ma¬ 
terial  possessing  unique  properties.  A  parametric  study  was  then  per¬ 
formed  by  varying  the  percentage  of  the  interface  modeled  as  perfectly 
weak. 

Again,  a  simplified  one  dimensional  mathematical  solution  and  both 
two  and  three  dimensional  finite  element  solutions  were  used  in  this  anal¬ 
ysis.  The  one  dimensional  model  is  a  modification  of  the  previous  model 
mentioned  in  the  last  section.  However,  an  additional  matrix  region  has 
been  added,  and  the  contact  between  it  and  the  fiber  is  assumed  to  be  per¬ 
fectly  weak(Figure  9).  Neglecting  shear  and  the  Poisson’s  effect  the  one 
dimensional  stress-strain  and  strain-displacement  equations  are  obtained 
as  before. 


Ef  =  !F=&  +  “f(T-Tp) 


(22) 


eii  = 


uir  uf  _  °i1 

b  El 


+  ai(T-Tp)  + 


£Pii 


(23) 
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Figure  9. 


Simplified  One  Dimensional  Model  With  Combined  Interphase 
Zone  /  Weak  Bond 
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(24) 


eMl.M  =  2£  +  aM(T-TP) 


Ei2=  ~  ^  4-  cti  (T  -  Tp)  +  epj2 


eM2=--2cUl1  =  ^  +  «m(T-Tp) 


=  ,~TT^=^  +  «m(T-Tp) 


EM2=Mbt^) 


The  three  equilibrium  equations  become 


of  a  +  Oi2  b  =  omi  d  +  o^  (a  +  b  -  d) 
Oij  (a  +  b  -  d)  =  OM2  (a  +  b  -  d) 

Omj  d  +  OM2  (a  +  b  -  d)  +  GM3  c  =  oa  (a  +  b  +  c) 


and  the  constraints  on  displacements 


ui2  “  ”F 

uMi  =  uM2  =  uM3 


The  same  procedure  is  followed  for  plastic  deformation  of  the  inter- 

orn  t»A«i  A«r  r»  r>  -i -r»  ♦v,  A  Alir  AAfi  A*.  U  ATITAtrAf  til  A  1»»A  aI/  U  A 
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tween  region  F  and  M 1  must  be  taken  into  account.  This  is  done  by  exam¬ 
ining  the  stress  in  region  M  1  after  each  load  step.  If  the  stress  is  com- 
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pressive,  then  uq  is  set  equal  to  up.  If  the  stress  in  region  M  \  is  tensile, 
then  this  stress  is  set  equal  to  zero  and  uq  becomes  the  unknown.  As  be¬ 
fore,  a  set  of  logic  sequences  are  followed  in  solving  the  equation  sets 
until  two  successive  solutions  are  identical.  The  Mathematica™  program 
incorporated  for  solving  this  combined  interphase  zone  /  weak  bond  sim¬ 
plified  one  dimensional  model  is  listed  in  appendix-C. 

Unlike  the  one  dimensional  model  where  only  the  percentage  of  inter¬ 
face  modeled  as  perfectly  weak  enters  into  the  calculations,  the  finite  ele¬ 
ment  solutions  also  are  affected  by  the  location.  A  choice  must  therefore 
be  made  as  to  where  to  place  the  perfectly  weak  bond  interface  in  the 
model.  For  both  the  two  and  three  dimensional  models  a  symmetric 
scheme  was  incorporated.  The  portion  of  the  interface  that  was  assumed 
to  behave  as  if  an  isotropic  interphase  zone  existed  between  the  fiber  and 
matrix  was  centered  about  a  45°  line  intersecting  the  center  and  upper 
right  corner  of  the  representative  volume  element  as  depicted  in 
Figure  10. 

The  grid  used  was  the  same  as  of  the  5%  of  fiber  radius  interphase 
zone  size  coarse  grid  model  shown  in  Figure  7.  Gap  elements  were  added 
symmetrically  as  the  size  of  the  interface  modeled  as  a  weak  bond  in¬ 
creased.  The  elements  outside  the  fiber  attached  to  the  gap  elements  were 
then  given  properties  identical  to  the  matrix. 

The  combined  interphase  zone  /  weak  bond  models  provide  a  good  ex¬ 
tension  on  the  effects  of  the  interface  conditions  on  the  overall  character¬ 
istics.  The  results  from  it  and  the  other  models  already  mentioned  are 
presented  in  the  following  chapters.  As  a  final  note,  no  mention  of  spe¬ 
cific  material  properties  given  to  a  particular  model  have  been  made  in 
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Figure  10.  Interphase  Zone  Location  Scheme  for  Combined  Finite  Ele¬ 
ment  Model 

the  previous  model  descriptions.  This  was  intentional.  The  following 
parametric  study  will  require  many  material  property  variations  to  be 
made.  Therefore,  it  was  thought  to  be  less  confusing  to  call  out  the  spe¬ 
cific  material  properties  used  in  a  particular  model  in  the  results  section. 
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IV.  Results  and  Discussion 


The  main  thrust  of  this  thesis  is  to  investigate  how  various  properties 
and  conditions  at  the  fiber/matrix  interface  in  a  high  temperature  metal 
matrix  composite  effect  its  overall  transverse  normal  properties.  The  fol¬ 
lowing  discussion  presents  the  results  of  such  a  parametric  study  complet¬ 
ed  by  using  the  techniques  revealed  in  the  previous  chapter.  All  thermo¬ 
mechanical  loading  sequences  consisted  of  cooldown  to  room  temperature 
from  a  specified  processing  temperature  followed  by  a  transverse  normal 
loading  at  constant  temperature  and  in  most  cases,  unloading.  The  trans¬ 
verse  stress-strain  response  during  the  mechanical  loading  and  unloading 
will  be  analyzed  in  detail.  First,  the  results  from  perfectly  strong  and 
perfectly  weak  interfaces  will  be  presented  and  compared  with  previous 
calculations  obtained  by  an  independent  researcher[2;  3].  Subsequently, 
the  results  from  models  containing  an  interphase  zone  between  fiber  and 
matrix  will  be  presented,  and  finally,  an  investigation  where  part  of  the 
interface  is  modeled  as  perfectly  weak  and  part  as  if  an  interphase  zone 
exists  will  be  put  forth. 

4.1  Strong  and  Weak  Interfacial  Bond  Results 

The  initial  portion  of  this  study  will  examine  the  behavior  of  a  finite 
element  micromechanics  model  where  the  interface  is  modeled  both  as 
perfectly  strong  and  perfectly  weak.  Since  previous  work  on  this  subject 
has  been  performed,  it  provides  an  excellent  means  for  validation  of  the 
calculations.  If  the  results  presented  in  this  section  show  good  agreement 
with  previous  work,  then  further  calculations  may  be  attempted  with  a 


high  degree  of  confidence.  Therefore,  much  of  the  material  properties 
and  other  parameters  will  be  matched  with  those  used  by  Nimmer[2;3]. 

4.1.1  Results  from  Perfectly  Strong  Interface 

As  described  in  section  3.3,  both  two  and  three  dimensional  finite  ele¬ 
ment  models  possessing  a  strong  interface  were  created.  A  35%  fiber  vol¬ 
ume  fraction  with  a  square  fiber  array  was  assumed,  and  the  calculations 
were  performed  using  the  material  properties  listed  in  Table  2.  One  of  the 
goals  during  the  strong  bond  calculations  was  to  determine  what  type  of 
two-dimensional  representation,  either  plane  stress  or  plane  strain,  more 
accurately  portrayed  a  generalized  plane  strain  analysis.  Therefore,  all 
three  methods  were  used  and  comparisons  made  between  the  plane  stress 
and  plane  strain  analyses  against  the  three-dimensional  generalized  plane 
strain  analysis.  The  loading  sequence  consisted  of  cooldown  from  a 
1000°C  zero  stress  state  to  room  temperature  at  23°C  followed  by  a  mono- 
tonically  increasing  transverse  normal  load.  The  Von  Mises  stresses  after 
cooldown  for  plane  stress  and  plane  strain  are  shown  in  Figure  1 1,  and  the 
generalized  plane  strain  model  contours  are  shown  in  Figure  12. 


Table  2.  Material  Properties  Used  for  Strong  Bond  Models  and  Weak 
Bond  RVE#1  Models 


Matrix 

Fiber 


E(GPa)  v  axlO^fC)-1 


— 

100.0 

0.3 

10.0 

400.0 

0.3 

5.0 

45 


ft  Figure  11.  Strong  Interface  Plane  Stress  and  Plane  Strain  Finite  Element 

Model  Residual  Stress  After  Cooldown 
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Figure  12.  Strong  Interface  Generalized  Plane  Strain  Finite  Element 
Model(3-D  Model)  Residual  Stress  After  Cooldown 


The  maximum  Von  Mises  stress  is  543  MPa  from  the  plane  stress  anal¬ 
ysis  and  745  MPa  from  the  plane  strain  analysis  while  the  generalized 
plane  strain  model  gives  794  MPa.  Also,  it  is  apparent  after  examination 
of  the  stress  contours  in  the  above  figures  that  from  a  stress  field  stand¬ 
point  the  plane  strain  analysis  more  accurately  approximates  the  three-di¬ 
mensional  generalized  plane  strain  model.  However,  since  the  goal  of  mi¬ 
cromechanics  is  to  approximate  the  overall  behavior,  the  load-deflection 
curve  of  the  entire  structure  is  of  primary  importance. 

The  overall  composite  horizontal  strain  is  the  displacement  of  the 
ripht-hand  face  referenced  to  the  after  cooldown  state  divided  hv  the  hori- 

LJ  '  v 

zontal  width,  and  the  overall  composite  stress  is  the  applied  load  on  the 
right-hand  face  divided  by  its  respective  area.  Thus,  the  overall  compos- 
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ite  transverse  normal  stress-strain  behavior  may  be  plotted  from  the  finite 
element  solutions.  This  is  done  in  Figure  13  for  plane  stress,  plane  strain, 
and  the  generalized  plane  strain  model,  and  the  same  procedure  will  be 
used  throughout  this  paper  for  plotting  all  transverse  response  curves. 

Surprisingly,  the  plane  stress  rather  than  the  plane  strain  solution 
agrees  more  closely  with  generalized  plane  strain.  As  expected,  plane 
stress  and  plane  strain  bounds  the  generalized  plane  strain  solution,  but 
the  plane  strain  solution  is  too  stiff.  Therefore,  when  considering  trans¬ 
verse  response  of  a  high  temperature  metal  matrix  composite,  a  plane 
stress  analysis  seems  to  be  the  best  two-dimensional  substitute  for  a 
three-dimensional  generalized  plane  strain  solution. 
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Figure  13.  Strong  Interface  Transverse  Response  for  Plane  Stress,  Plane 
Strain,  and  Generalized  Plane  Strain  Finite  Element  Solutions 
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A  possible  explanation  for  this  apparent  anomaly  lies  with  the  thermal 
strains.  The  thermal  strains,  which  are  not  shown  in  the  figure  since  the 
zero  strain  state  is  taken  to  be  the  position  of  the  right-hand  face  after 
cooldown,  are  on  the  order  of  10-3  while  the  mechanical  strains  are  on  the 
order  of  10-4.  Thus,  the  thermal  strains  are  an  order  of  magnitude  greater 
than  the  mechanical  strains,  and  since  thermal  strains  are  triaxial  in  na¬ 
ture,  a  significant  amount  of  contraction  will  occur  in  the  out-of-plane  di¬ 
rection.  Hence,  the  type  of  out-of-plane  constraint  will  have  a  pro¬ 
nounced  effect  on  the  behavior.  A  pictorial  explanation  is  given  in 
Figure  14.  A  generalized  plane  strain  analysis  allows  uniform  out-of- 


Figure  14.  Schematic  of  Out-of-Plane  Constraints  Showing  Difficulty  of 
a  Plane  Strain  Solution  in  the  Presence  of  Large  Thermal 
Strains 
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plane  longitudinal  contraction  during  cooldown  while  the  plane  stress 
analysis  allows  free  longitudinal  contraction.  On  the  other  hand,  if  a 
plane  strain  solution  is  used,  no  thermal  contraction  is  allowed  in  the  lon¬ 
gitudinal  direction,  but  this  contraction  is  large  when  compared  to  the  me¬ 
chanical  displacements.  Therefore,  the  plane  strain  constraint  will  result 
in  unduly  large  stresses  in  the  longitudinal  direction.  The  Poisson’s  ef¬ 
fect  of  these  longitudinal  stresses  will  produce  a  stiffer  structure  in  the 
lateral  or  transverse  direction. 

This  last  observation  is  important  from  the  computational  aspect  since 
three-dimensional  nonlinear  finite  element  solutions  are  costly  in  CPU 
time.  For  the  strong  interface  models  where  very  little  nonlinearity  oc¬ 
curs  it  was  possible  to  use  a  linear  analysis,  but  for  all  other  models  a 
nonlinear  analysis  was  necessary.  Therefore,  appropriate  two-dimensional 
solutions  must  be  sought.  Nimmer  used  a  plane  stress  type  solution  while 
performing  his  two-dimensional  finite  element  analysis,  and  therefore,  if 
appropriate  comparisons  are  to  be  made  with  his  results,  then  plane  stress 
solutions  must  be  used.  In  addition,  a  programming  error  presently  exists 
in  MSC/NASTRAN  for  solutions  involving  a  plane  strain  nonlinear  mate¬ 
rial  analysis  in  the  presence  of  thermal  loads.  Hence,  for  these  reasons, 
in  subsequent  sections  all  two-dimensional  models  will  use  a  plane  stress 
type  solution. 

4.1.2  Results  from  Perfectly  Weak  Interface 

The  finite  element  solutions  of  a  strong  interface  provided  a  good 
comparison  to  the  weak  interface  models  since  before  fiber/matrix  separa¬ 
tion  occurs  under  transverse  loading  the  weak  and  strong  bond  models 
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should  behave  in  the  same  manner.  Since  nonlinear  solutions  must  be 
used,  the  weak  interface  greatly  increases  the  complexity  of  the  problem. 
Therefore,  to  be  assured  that  the  results  were  reasonable  several  calcula¬ 
tions  were  performed  using  a  weak  bond. 

Initially,  a  35%  fiber  volume  fraction  representative  volume  element 
(RVE#1)  similar  to  that  used  for  a  strong  bond  was  incorporated.  The 
loading  sequence  for  the  RVE#1  weak  bond  models  was  the  same  as  for 
the  strong  bond  models.  The  overall  composite  transverse  stress-strain  re¬ 
sponse  after  cooldown  for  a  monotonically  increasing  load  using  a  plane 
stress  analysis  is  given  in  Figure  15  and  compared  with  a  strong  bond. 

The  material  properties  were  the  same  as  listed  in  Table  2  on  page  45,  and 
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Figure  15.  Weak  Interface  Transverse  Response  of  RVE#1  Assuming 
Linear  Elastic  Matrix  and  Using  a  Plane  Stress  Finite  Ele¬ 
ment  Solution 
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the  matrix  was  assumed  to  be  linearly  elastic.  Also,  a  Coulomb  friction 
coefficient  of  0.0  was  input  for  the  gap  elements  along  the  interface.  The 
classical  bilinear  shape  of  the  weak  interface  curve  is  demonstrated,  and 
fiber/matrix  sep.  ration  occurs  at  approximately  150  MPa. 

Notwithstanding,  the  initial  slope  of  the  weak  bond  curve  does  not 
match  that  of  the  strong  bond  as  expected.  An  investigation  into  this  dif¬ 
ference  requires  an  examination  of  the  displaced  geometry  (Figure  16). 

The  deformed  structure  is  given  for  a  transverse  load  of  140  MPa  and 
210  MPa.  The  140  MPa  load  point  occurs  before  separation,  but  upon  fur¬ 
ther  examination  of  the  displaced  geometry  it  is  apparent  that  there  is  an¬ 
other  mechanism  acting  besides  just  fiber/matrix  separation.  A  weak  in¬ 
terface  with  no  friction  will  also  experience  slip  along  the  interface. 
Therefore,  even  though  the  fiber  and  matrix  are  still  in  contact  at 
140  MPa,  the  structure  has  been  considerably  weakened  as  opposed  to  the 
strong  interface  due  to  the  fiber/matrix  slip.  Interfacial  slip  was  found  to 
occur  at  all  load  levels,  and  hence,  the  decreased  initial  slope  in  the 
stress-strain  curve  resulted. 

As  a  side  note,  the  plots  in  this  paper  that  show  deformed  geometry 
are  presented  with  deformations  grossly  exaggerated,  usually  by  multipli¬ 
cation  factors  above  100.  Many  of  the  stress  contour  plots  are  shown  with 
their  deformed  geometry. 

The  addition  of  Coulomb  friction  at  the  interface  has  a  pronounced  ef¬ 
fect  on  the  results  since  it  tends  to  arrest  the  fiber/matrix  slip  (Figure  17). 
As  the  friction  coefficient  is  increased,  the  initial  modulus  before 
fiber/matrix  separation  also  increases  until  it  becomes  equal  to  that  of  a 
strong  bond.  However,  the  addition  of  friction  has  no  effect  on  the  sec- 
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0.777X10-4  mm 
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of  210  MPa 


Note:  Fiber  Radius  =  .07  mm 
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Figure  16.  Weak  Interface  Deformed  Geometry  of  RVE#1  at  140  MPa 
and  210  MPa. 


53 


Strain 

Figure  17.  Weak  Interface  Transverse  Response  of  RVE#1  With  and 

Without  Friction  at  the  Interface  and  Assuming  a  Linear  Elas¬ 
tic  Matrix 


ondary  slope  that  occurs  after  separation.  Since  a  large  portion  of  the 
fiber  and  matrix  surfaces  are  no  longer  in  contact,  the  frictional  constraint 
vanishes. 

The  final  calculation  using  a  monotonically  increasing  load  incorpo¬ 
rated  an  elastic-perfectly  plastic  matrix  with  a  yield  stress  of  800  MPa 
and  a  friction  coefficient  of  0.0  at  the  interface(Figure  18).  All  other  ma¬ 
terial  properties  remained  the  same  as  in  Table  2.  The  results  are  plotted 
alongside  that  obtained  by  Nimmer  using  the  finite  element  program 
ADINA  for  comparison^].  As  shown,  the  agreement  between  the  calcula¬ 
tions  is  excellent.  Especially  when  one  considers  that  the  grids  were  cre- 
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Figure  18.  Weak  Interface  Transverse  Response  of  RVE#1  Assuming 

Elastic-Plastic  Matrix  With  a  Yield  Stress  of  800  MPa  and  No 
Interfacial  Friction  as  Compared  to  Nimmer’s  Results  [2] 

ated  independently.  Fiber/matrix  separation  again  occurs  at  approximate¬ 
ly  150  MPa,  but  the  transverse  response  after  separation  is  not  linear  since 
the  matrix  exhibits  yielding. 

A  more  clear  understanding  of  the  effects  of  the  monotonically  in¬ 
creasing  load  may  be  obtained  by  examining  the  history  of  the  stress  field 
during  loading.  The  Von  Mises  stress  contours  after  cooldown  and  at  the 
onset  of  fiber/matrix  separation  are  presented  in  Figure  19.  The  stress  for 
transverse  loads  of  210  MPa  and  260  MPa  are  then  shown  in  Figure  20. 
The  maximum  Von  Mises  stress  after  cooldown  is  621  MPa.  Therefore, 
for  a  matrix  yield  stress  of  800  MPa  the  residual  stresses  during  cooldown 
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Figure  19.  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#1  With 
Elastic-Plastic  Matrix  and  Frictionless  Interface  at  Cooldown 
and  Onset  of  Fiber/Matrix  Separation 
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Figure  20.  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#1  With 
Elastic-Plastic  Matrix  and  Frictionless  Interface  at  210  MPa 
and  260  MPa  Transverse  Load 
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are  not  sufficient  to  cause  matrix  plastic  flow.  Also,  when  compared  to 
Figure  1 1  on  page  46  the  stress  field  after  cooldown  for  a  weak  and  strong 
bond  are  not  identical  as  one  might  at  first  imagine.  The  weak  interface 
with  no  friction  is  unable  to  transmit  shear  across  its  face  since  friction  is 
absent,  so  slip  occurs  and  causes  what  seems  to  be  an  apparent  discontinu¬ 
ity  at  the  interface  in  the  stress  contours.  Large  compressive  stresses  nor¬ 
mal  to  the  interface  exist  after  cooldown.  These  compressive  stresses  are 
maximum  along  the  fixed  edges  as  shown  in  Figure  19. 

Interfacial  compressive  stresses  must  be  overcome  before  the  fiber 
and  matrix  will  separate.  This  separation  occurs  at  150  MPa,  and  as 
shown  in  the  figure  a  very  small  portion  of  the  matrix  has  exceeded  the 
yield  stress  at  this  point,  indicating  that  for  the  model  considered,  matrix 
yield  occurs  simultaneously  with  fiber/matrix  separation.  This  was  mere¬ 
ly  a  coincident.  In  general,  fiber/matrix  separation  and  matrix  yield  are 
not  required  to  occur  at  the  same  load  level. 

As  the  model  is  further  loaded  past  the  separation  point,  plasticity  in 
the  matrix  grows.  Also,  as  more  of  the  matrix  separates  from  the  fiber  the 
overall  stress  in  the  fiber  begins  to  decrease  (Figure  20). 

Weak  interface  calculations  involving  a  loading  and  unloading  se¬ 
quence  used  the  32%  fiber  volume  fraction,  RVE#2,  finite  element  grid. 
Both  plane  stress  and  three-dimensional  generalized  plane  strain  solutions 
were  performed.  As  mentioned  in  chapter  3,  this  finite  element  model 
was  created  to  match  published  experimental  data[3].  Therefore,  specific 
material  properties  associated  with  the  published  material  were  used  and 
are  given  in  Table  3.  A  coefficient  of  friction  at  the  interface  of  0.30  was 
chosen  to  match  Nimmer’s  results,  and  the  loading  sequence  consisted  of 
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Table  3.  Material  Properties  Used  for  Weak  Bond  RVE#2  Models 


E(GPa) 

V 

a  xio-ficc)-* 

Gy$  (MPa) 

H'  (MPa) 

Matrix 

113.7 

0.3 

9.44 

900.0 

4600 

Fiber 

414.0 

0.3 

4.86 

- 

- 

cooldown  from  a  900°C  zero  stress  state  to  21°C  followed  by  a  transverse 
applied  load  at  constant  temperature  to  425  MPa  and  unload[3]. 

The  transverse  response  after  cooldown  for  both  plane  stress  and  gen¬ 
eralized  plane  strain  is  presented  in  Figure  21.  Results  from  Nimmer’s 
generalized  plane  strain  analysis  are  plotted  alongside  for  comparison. 

The  difference  between  the  plane  stress  and  generalized  plane  strain  re¬ 
sults  are  significant.  The  point  of  fiber/matrix  separation  is  slightly  in¬ 
creased  and  the  overall  response  is  considerably  stiffer  for  the  generalized 
plane  strain  model.  However,  the  initial  moduli  are  approximately  the 
same,  so  the  deviation  occurs  after  the  separation  point.  This  deviation 
between  the  generalized  plane  strain  solution  and  plane  stress  is  linked  to 
the  high  triaxial  state  of  stress  that  exists  for  generalized  plane  strain  due 
to  the  extra  constraint  in  the  out-of-plane  direction.  High  triaxial  stress 
fields  may  still  possess  low  octahedral  shear  stress  which  may  not  be  suf¬ 
ficient  to  cause  yielding.  Hence,  if  a  triaxial  stress  field  is  the  cause  of 
the  large  difference  between  the  generalized  plane  strain  and  plane  stress 
solutions,  then  a  significant  reduction  in  matrix  yield  should  be  visible  in 
the  stress  contours  (Figure  22).  Such  is  found  to  be  the  case  as  the 
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Figure  21.  Weak  Interface  Transverse  Response  for  RVE#2  Using  Plane 
Stress  and  Generalized  Plane  Strain  Solutions  With  Compari¬ 
son  to  Nimmer’s  Results  [3] 


900  MPa  contour  line  which  represents  the  extent  of  matrix  yield  is  great¬ 
ly  reduced  in  area  for  the  generalized  plane  strain  solution. 

Nevertheless,  as  pointed  out  early  on  in  Figure  2  on  page  16  the  weak 
bond  generalized  plane  strain  solutions  do  not  always  match  experiments. 
Hence,  it  is  possible  that  additional  constraints  exist  in  the  actual  speci¬ 
men  that  are  not  modeled  with  a  completely  weak  interfacial  bond.  A 
partially  weak  bond  is  more  physically  realistic,  but  the  finite  element 
weak  bond  model  presented  in  the  literature [2;  3]  and  also  up  to  this  point 
in  the  present  study  offers  little  flexibility  for  incorporating  additional 
constraints;  interfacial  friction  is  the  only  variable.  An  increase  in  the  in- 
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Figure  22.  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#2  at 

Maximum  Load  for  Both  Plane  Stress  and  Generalized  Plane 
Strain 
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terfacial  coefficient  of  friction  was  previously  found  to  increase  the  mod¬ 
ulus  before  separation  but  had  little  or  no  effect  after  separation. 

However,  the  problem  at  hand  requires  a  constraint  that  exists  both  before 
and  after  separation. 

The  zone  between  fiber  and  matrix  is  very  complex.  During  high  tem¬ 
perature  processing  a  reaction  zone  is  formed  where  the  fiber  material  has 
penetrated  into  and  reacted  with  the  matrix.  Many  flaws  exist  in  this  zone 
and  its  failure  occurs  early  during  initial  loading  [l].  After  failure  of  the 
fiber/matrix  bond  a  fracture  surface  exists.  Most  fracture  surfaces  are  not 
smooth  but  rather  jagged  in  nature,  especially  at  the  microscopic  level. 
Therefore,  the  surface  will  possess  a  high  resistance  to  slip  between  the 
two  faces,  and  as  long  as  the  separation  between  the  two  fracture  surfaces 
is  small,  this  resistance  to  slip  will  persist. 

Hence,  a  reasonable  constraint  to  add  to  the  perfectly  weak  bond 
model  is  one  that  prevents  slip  between  fiber  and  matrix  both  before  and 
after  separation.  Such  a  constraint  was  added  to  the  plane  stress  and  gen¬ 
eralized  plane  strain  models  to  determine  its  effect,  figure  23  presents 
the  plane  stress  solution.  The  added  constraint  of  preventing  slip  at  the 
interface  has  a  pronounced  effect  on  the  result,  greatly  stiffening  the 
structure.  In  fact,  it  increases  the  overall  stiffness  to  the  point  that  no 
plastic  deformation  occurs  in  the  matrix  making  the  loading  and  unloading 
elastic.  The  initial  modulus  remains  the  same,  and  fiber/matrix  separation 
occurs  at  approximately  the  same  point,  but  the  slope  after  separation  is 
greatly  increased. 

The  generalized  plane  strain  solution  is  shown  in  Figure  24  with  the 
experimental  results  from  Figure  2  on  page  16  plotted  alongside  for  com- 
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Figure  23.  Weak  Interface  Transverse  Response  for  RVE#2  Plane  Stress 
Solution  With  and  Without  Fiber/Matrix  Slip  at  the  Interface 


parison.  In  addition,  the  stress  contours  at  maximum  load  are  given  in 
Figure  25.  The  results  again  indicate  a  substantial  effect  from  the  addi¬ 
tional  constraint.  In  the  same  way  as  in  the  plane  stress  solution  the 
structure  has  been  greatly  stiffened.  In  fact,  the  Von  Mises  stress  of  the 
matrix  never  exceeds  its  yield  stress  of  900  MPa.  Instead,  the  overall 
transverse  stress-strain  response  becomes  a  nonlinear  elastic  curve. 
Fiber/matrix  separation  still  occurs,  but  due  to  the  extra  constraint  that 
prevents  slip,  the  overall  transverse  response  is  dramatically  changed,  and 
ne  separation  distance  at  maximum  load  is  greatly  reduced.  Also,  a  very 
encouraging  development  is  that  the  with  and  without  slip  constraint  cal¬ 
culations  form  a  bound  for  the  experimental  data. 
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Figure  24.  Weak  Interface  Transverse  Response  for  RVE#2  General¬ 
ized  Plane  Strain  Solution  With  and  Without  Fiber/Matrix 
Slip  and  Comparison  to  Nimmer’s  Experiment^] 


As  a  reminder,  the  experimental  data  is  from  a  load,  unload,  and  re¬ 
load  sequence.  The  difference  between  the  initial  load  and  reload  curves 
in  the  experimental  data  are  from  the  finite  fiber/matrix  bond  strength 
which  fails  during  initial  loading.  The  experiment  exhibits  energy  dissi¬ 
pation  during  the  cycle,  but  very  little  permanent  strain.  Therefore,  it  is 
reasonable  to  assume  that  much  of  the  energy  dissipation  is  caused  by 
something  other  than  plastic  deformation.  During  the  opening  and  closing 
of  jagged  fracture  surfaces  frictional  losses  will  occur  that  otherwise 
would  not  if  the  two  surfaces  were  smooth.  Hence,  the  motivation  for  the 
additional  constraint  also  serves  as  an  explanation  for  the  energy  dissipa¬ 
tion  observed  in  experiment.  It  was  not  possible  to  account  for  this  fric- 
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Figure  25.  Weak  Interface  Von  Mises  Stress  Contours  for  RVE#2  at 

Maximum  Load  Using  Generalized  Plane  Strain  Solution  and 
Preventing  Fiber/Matrix  Slip 

tional  energy  loss  in  the  calculations  as  that  would  almost  certainly  in¬ 
volve  the  development  of  a  new  type  of  element  in  MSC/NASTRAN. 

The  assumption  of  a  jagged  fiber/matrix  fracture  surface  provides  a 
possible  explanation  for  the  phenomena  observed.  However,  as  demon¬ 
strated,  incorporating  its  effects  in  an  analysis  is  difficult.  In  spite  of  the 
jagged  surface  a  certain  amount  of  fiber/matrix  slip  will  still  be  allowed, 
and  frictional  losses  during  opening  and  closing  make  the  problem  even 
more  formidable.  At  present,  it  is  impossible  to  include  such  effects  with 
the  finite  element  tools  available.  Therefore,  an  attempt  is  made  to  model 
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the  problem  in  a  simplified  manner  by  including  a  third  zone  between 
fiber  and  matrix  possessing  its  own  unique  properties.  This  interphase 
zone  should  provide  for  energy  dissipation  and  for  additional  constraints 
at  the  interface  not  modeled  by  a  perfectly  weak  bond. 

4.2  Results  from  Models  Possessing  an  Interphase  Zone 

Many  mechanisms  are  acting  at  the  fiber/matrix  interface.  Imperfec¬ 
tions  exist  that  have  a  pronounced  effect  on  the  composite’s  overall  char¬ 
acteristics.  The  challenge  is  to  model  these  mechanisms  in  the  simplest 
way.  As  demonstrated  in  the  previous  section  assuming  the  interface  be¬ 
haves  as  if  it  is  perfectly  weak  exhibits  modest  success,  but  due  to  inade¬ 
quate  constraints  at  the  interface  and  because  there  is  no  means  of  model¬ 
ing  energy  dissipation  such  an  assumption  meets  with  difficulty.  In  this 
section  the  attempt  is  made  to  model  the  interface  as  a  third  zone  possess¬ 
ing  unique  properties.  Previous  studies  have  used  such  an  approach,  but 
they  only  examined  the  elastic  realm[17-19].  This  study  will  emphasize 
the  plastic  properties.  The  assumption  being  that  all  the  nonlinearities 
existent  at  the  interface  may  be  adequately  modeled  by  an  isotropic  elas¬ 
tic-plastic  interphase  zone  between  the  fiber  and  matrix. 

4.2.1  One-Dimensional  Results  With  Interphase  Zone 

The  simplified  one-dimensional  model  discussed  in  section  3.4  was 
incorporated  for  analyzing  parameter  variations  associated  with  an  inter-  - 
phase  zone.  The  advantage  of  the  one-di  nensional  solutions  over  the  fi¬ 
nite  element  calculations  rests  with  the  speed  at  which  answers  could  be 


obtained.  Solutions  that  took  on  the  average  ninety  seconds  with  the  one- 
dimensional  model  would  require  several  hours  for  a  finite  element  solu¬ 
tion.  In  addition,  the  finite  element  solutions  would  generally  require  the 
creation  of  a  new  grid  which  in  itself  might  take  several  hours.  There¬ 
fore,  the  mode  of  analysis  presented  here  proved  to  be  very  helpful  in  this 
study. 

Due  to  the  fact  that  several  parameter  variations  were  to  be  examined 
a  suitable  baseline  was  chosen  from  which  a  single  parameter  could  be 
varied  while  holding  the  others  fixed.  This  baseline  set  is  listed  in 
Table  4  and  was  chosen  to  match  the  material  properties  from  some  of  the 
weak  bond  calculations  of  the  previous  section.  All  parameters  associated 
with  the  interphase  zone  were  varied  to  determine  their  overall  effect. 
These  include  interphase  zone  yield  stress,  strain  hardening,  Young’s 

Table  4.  Baseline  Parameters  for  One-Dimensional  Model  With  Inter¬ 
phase 


Gys] 

-  125  MPa 

vF 

=  0.35 

H' 

-  0  MPa 

vb 

=  0.07  mm 

(see  Figure  6) 

TP 

=  iooo°c 

=  5%  of  fiber 

T 

=  25°C 

EF 

-  400  GPa 

CCF 

=  5  x  10-6(°C)-1 

El 

=  100  GPa 

aI 

=  10  x  10-6  (°C)-1 

Em 

=  100  GPa 

=  10  x  10-6  (°C)*1 
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modulus,  thermal  coefficient  of  expansion,  and  thickness.  Each  solution 
consisted  of  transverse  normal  loading  to  27.'  MPa  followed  by  a  complete 
unloading. 

Some  of  the  results  are  compared  to  a  one-dimensional  model  with  a 
strong  and  weak  bond[2].  This  is  the  case  in  Figure  26  where  interphase 
zone  yield  stress  is  varied  from  200  MPa  to  75  MPa.  As  shown,  the  initial 
slopes  of  all  curves  match  the  strong  interfacial  bond  slope  exactly.  This 
is  as  expected  since  the  initial  compressive  thermal  stresses  at  the  inter¬ 
face  must  be  overcome  before  nonlinearities  may  occur.  The  weak  bond 
presents  a  bilinear  behavior  where  the  “kink”  in  the  curve  results  from 
fiber/matrix  separation.  The  model  containing  an  elastic-plastic  inter- 
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Figure  26.  One-Dimensional  So  Jtions  With  Interphase  Zone  for 

Interphase  Yield  s' tier- cs  from  75  MPa  to  200  MPa  as  Com¬ 
pared  to  Nijnr-er’s  Stionv,  and  Weak  Bond  Model[2] 
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phase  zone  must  further  be  loaded  before  nonlinearity  develops  because 
the  interphase  yield  stress  must  also  be  overcome.  Similarly,  unloading  is 
characterized  by  elastic  followed  by  plastic  deformation  which  involved 
compressive  yielding  of  the  interphase.  The  plastic  deformation  during 
unloading  is  a  result  of  the  thermal  residual  stresses  and  elastic  energy  in 
the  matrix  overcoming  the  interphase  zone  yield  stress.  As  shown,  for  the 
given  loading  condition  a  yield  stress  of  200  MPa  does  not  experience 
yielding,  and  hence  a  linear  elastic  curve  that  matches  the  strong  bond  re¬ 
sults. 

At  first  one  might  reason  that  a  yield  stress  of  zero  should  match  the 
weak  bond  model,  but  in  actuality  a  zero  yield  stress  results  in  a  complete 
relaxation  of  all  thermal  residual  stresses  as  the  interphase  zone  deforms 
during  cooldown.  The  compressive  stress  at  the  interphase  after  cooldown 
is  130  MPa.  Therefore,  any  yield  stress  below  this  will  result  in  some  re¬ 
laxation  of  the  thermal  stresses.  This  was  observed  in  the  solutions  and 
c.n  be  visualized  in  Figure  26.  For  yield  stresses  above  130  MPa  the 
point  where  nonlinearity  first  occurs  is  controlled  by  the  applied  load 
overcoming  an  initial  compressive  s+ress  of  130  MPa  at  the  interphase  and 
then  exceeding  the  tensile  yield  stress.  For  yield  stresses  below  130  MPa 
the  initial  compressive  stress  that  the  applied  load  must  overcome  is  equal 
to  the  yield  stress.  Hence,  the  point  of  nonlinearity  will  occur  sooner. 

This  can  be  seen  by  observing  that  all  curves  1th  yield  stresses  below 
130  MPa  demonstrate  a  greater  difference  between  points  of  nonlinearity 
for  two  successive  yield  stresses  than  can  be  seen  in  the  curves  with  yield 
stiesses  above  130  MPa.  Therefore,  rather  than  zero  yield  stress  matching 
the  weak  bond,  a  yield  stress  of  75  MPa  closely  approximates  the  behavior 


of  a  weak  bond  during  loading.  Thus,  an  incorrect  stress  held  calculation 
will  result  due  to  the  relaxation  of  the  thermal  residual  stresses. 

In  addition,  the  simplified  one-dimensional  solutions  demonstrate  that 
a  significant  amount  of  energy  dissipation  occurs  when  modeling  the 
fiber/matrix  interface  as  a  third  zone  possessing  relatively  low  yield 
stress.  The  area  traced  out  by  each  loading  and  unloading  sequence  de¬ 
picts  this  energy  dissipation.  Also,  for  yield  stresses  of  125  MPa  and 
below  the  model  returns  close  to  its  approximate  initial  strain  state  after 
unloading.  This  results  from  the  residual  thermal  compressive  stresses 
being  equal  to  the  yield  stress,  and  hence  they  are  capable  of  plastically 
deforming  the  interphase  zone  back  to  its  original  state.  Therefore,  this 
approach’s  ability  to  model  energy  dissipation  and  allow  for  additional  in¬ 
terface  constraints  that  may  be  varied  to  approximate  any  level  of  bonding 
from  strong  to  weak  holds  some  promise  over  the  weak  bond  model. 

Varying  the  plastic  parameters(i.e.  yield  stress  and  strain  hardening) 
was  found  to  effect  the  overall  transverse  behavior  the  most.  The  effect 
of  strain  hardening  on  the  solution  is  shown  in  Figure  27,  and  it  indicates 
that  the  main  effect  is  on  the  secondary  slope  during  loading.  As  a  result, 
when  the  interphase  strain  hardening  parameter  is  increased  the  energy 
dissipated  during  loading  and  unloading  is  decreased.  Also,  a  mirror 
image  of  these  slopes  about  the  H'=  0  MPa  secondary  slope  line  occurs 
when  plastic  deformation  becomes  present  upon  unloading.  This  is  a  re¬ 
sult  of  the  isotropic  strain  hardening  assumption.  In  fact,  for  the 
H'=  2000  MPa  curve  the  interphase  has  been  sufficiently  hardened  during 
loading  that  the  residual  thermal  stresses  are  unable  to  induce  plastic  de¬ 
formation  when  the  structure  is  unloaded. 
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Figure  27.  One-Dimensional  Solutions  of  Interphase  Zone  Strain  Hard¬ 
ening  Parameter  Variation 

The  three  other  parameters  varied  were  the  interphase  zone  size, 
Young’s  modulus,  and  thermal  coefficient  of  expansion.  The  effect  of 
zone  size  is  depicted  in  Figure  28.  As  shown,  it  has  only  modest  impact 
on  the  results.  A  slight  change  in  the  secondary  slope  occurs  reducing  the 
amount  of  energy  dissipated  as  the  zone  size  is  decreased.  This  is  as  ex¬ 
pected  since  the  matrix  is  stronger  than  the  interphase  7 one.  Therefore,  as 
the  interphase  zone  size  is  decreased  and  replaced  by  matrix  material  the 
overall  structure  will  become  stronger. 

Young’s  modulus  was  also  found  to  have  only  a  slight  effect  on  the 
overall  stress-strain  response  during  loading  and  unloading(Figure  29). 
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Figure  75.  One-Dimensional  Solutions  of  Interphase  Zone  Size  Variation 


Since  the  interphase  zone  is  small  in  comparison  to  the  fiber  and  matrix, 
its  elastic  characteristics  exert  little  influence  on  the  entire  structure.  In 
the  figure  below  the  interphase  modulus  is  first  set  equal  to  that  of  the 
matrix(100  GPa)  and  then  the  fiber(400  GPa).  Only  a  slight  change  in  the 
elastic  slope  of  the  structure  results. 

Finally,  the  thermal  coefficient  of  expansion  of  the  interphase  is  var¬ 
ied  and  shown  in  Figure  30.  As  with  the  modulus,  the  expansion  coeffi¬ 
cient  is  varied  from  the  the  matrix’s  value  to  the  fiber’s.  The  plot  shown 
deviates  from  the  baseline  parameter  set  by  taking  the  interphase  yield 
stress  to  be  150  MPa.  This  was  necessary  since  the  effect  of  changing  the 
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Figure  29.  One-Dimensional  Solutions  of  Interphase  Young’s  Modulus 
Variation 

expansion  coefficient  is  reflected  in  the  compressive  residual  stresses,  and 
in  order  to  prevent  a  relaxation  of  these  stresses  a  yield  stress  was  chosen 
so  that  yielding  would  not  occur  during  cooldown.  With  transverse  load¬ 
ing  as  the  expansion  coefficient  is  increased  the  stress  at  which  nonlinear¬ 
ity  begins  decreases.  This  results  from  the  interphase  zone  with  a  higher 
expansion  coefficient  contracting  more  during  cooldown  and  hence  the 
matrix  experiences  less  resistance  when  it  contracts  thereby  lowering  the 
compressive  residual  stresses.  However,  this  effect  is  small. 

In  summary,  the  one-dimensional  solutions  have  provided  a  good* in¬ 
troduction  into  how  the  properties  of  the  interphase  zone  effect  the  over- 
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Figure  30.  One-Dimensional  Solutions  of  Interphase  Zone  Thermal  Ex¬ 
pansion  Coefficient  Variation 

all  characteristics  of  the  composite.  Yield  stress  and  strain  hardening 
were  found  to  have  the  most  profound  effect.  Increasing  the  yield  stress 
was  found  to  increase  the  point  of  nonlinearity  (Figure  26),  and  increasing 
the  strain  hardening  parameter  was  found  to  increase  the  secondary  slope 
(Figure  27).  Significant  energy 'dissipation  occurs  during  a  loading  and 
unloading  sequence,  and  for  sufficiently  low  interphase  zone  yield  stress 
the  overall  structure  returns  to  its  approximate  original  strain  state  after 
unloading.  A  more  thorough  investigation  using  finite  element  analysis  is 
presented  in  the  following  section,  and  some  comparisons  will  be  made 
between  it  and  the  one-dimensional  solutions. 
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L2.2  Finite  Element  Solutions  With  Interphase  Zone 

An  examination  into  the  effects  of  an  interphase  zone  between  fiber 
and  matrix  using  finite  element  analysis  utilized  the  finite  element  grids 
presented  in  section  3.4  for  a  35%  fiber  volume  fraction  and  a  fiber  radius 
of  0.07  mm.  The  bulk  of  the  calculations  consisted  of  plane  stress  type 
solutions,  but  a  few  three-dimensional  generalized  plane  strain  solutions 
were  also  performed. 

Just  as  for  the  one-dimensional  solutions,  a  baseline  parameter  set  was 
chosen  from  which  a  single  parameter  could  be  varied  while  holding  the 
others  fixed.  This  baseline  set  is  given  in  Table  5,  and  except  for  the 
added  matrix  plastic  properties,  it  is  the  same  as  was  used  in  the  one-di¬ 
mensional  solutions.  The  loading  sequence  is  also  the  same,  cooldown 
from  a  1000°C  processing  temperature  to  room  temperature  at  25°C  fol¬ 
lowed  by  transverse  loading  at  constant  temperature  to  275  MPa  and  then 
unload. 


Tabic  5.  Baseline  Parameters  for  Finite  Element  Models  With  Inter¬ 
phase  Zone 


E(GPa; 

V 

a  xio-6^)- 

l  Oys  (MPa) 

H'  (MPa) 

Matrix 

100.0 

0.3 

10.0 

800.0 

0. 

Interphase 

100.0 

0.3 

10.0 

125.0 

0. 

Fiber 

400.0 

0.3 

5.0 

- 

- 

Interphase  Zone  Size(thickness)  =  5%  of  Fiber  Radius 


The  plane  stress  solutions  will  be  presented  first.  The  effect  of  vary¬ 
ing  the  yield  stress  is  shown  in  Figure  31  along  with  the  solutions  for  a 
perfectly  strong  bond  and  a  perfectly  weak  bond.  The  behavior  is  very 
similar  to  that  seen  in  the  one-dimensional  solutions.  The  initial  modulus 
for  all  curves  matches  that  of  a  strong  bond,  and  the  main  effect  of  in¬ 
creasing  the  interphase  zone  yield  stress  is  increasing  the  stress  at  which 
nonlinearity  first  occurs.  An  interphase  yield  stress  of  200  MPa  experi¬ 
ences  very  little  yielding,  and  hence  this  curve  does  not  deviate  much 
from  the  strong  bond  solution. 

Relaxation  of  the  thermal  residual  stresses  during  cooldown  occurred 
in  the  one-dimensional  solutions  for  interphase  yield  stresses  below 
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Figure  31.  Finite  Element  Plane  Stress  Solutions  of  Interphase  Zone 

Yield  Stress  Variation  as  Compared  to  Perfectly  Strong  and 
Perfectly  Weak  Bonds 
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130  MPa.  A  more  thorough  examination  into  this  problem  by  using  the 
stress  fields  from  the  finite  element  solutions  may  provide  greater  insight. 
The  Von  Mises  stress  contours  after  cooldown  for  a  strong  bond  and  an 
interphase  yield  stress  of  175  MPa  are  given  in  Figure  32.  The  stress  con¬ 
tours  for  interphase  yield  stresses  of  125  MPa  and  75  MPa  are  given  in 
Figure  33.  One  can  see  that  a  relaxation  of  the  stresses  from  a  perfectly 
strong  bond  stress  field  occurs  to  a  certain  extent  in  all  three  solutions 
shown  containing  an  interphase  zone.  Of  course,  the  most  dramatic  effect 
occurs  with  the  75  MPa  interphase  yield  stress  solution.  The  maximum 
Von  Mises  stress  after  cooldown  for  the  strong  interface  is  545  MPa.  If 
an  interphase  zone  with  a  yield  stress  of  75  MPa  is  used,  this  maximum 
stress  is  reduced  to  345  MPa.  Therefore,  just  as  in  the  one-dimensional 
solutions  a  relaxation  of  the  thermal  stresses  also  occurs  in  the  finite  ele¬ 
ment  solutions,  but  in  a  more  gradual  fashion.  Portions  of  the  interphase 
zone  in  the  finite  element  solutions  are  at  a  relatively  high  state  of  stress 
after  cooldown.  All  yield  stresses  examined  were  below  this  maximum 
stress,  and  therefore  yielding  would  occur  in  these  regions.  On  the  other 
hand,  the  one-dimensional  solutions  required  the  stress  throughout  the  in¬ 
terphase  to  be  constant,  and  hence  the  entire  interphase  zone  was  either 
above  or  below  the  yield  stress. 

A  major  difference  between  the  one-dimensional  and  finite  element 
solutions  was  that  for  all  interphase  zone  yield  stresses  using  the  finite  el¬ 
ement  analysis  the  composite  experienced  permanent  strain  after  unload¬ 
ing.  However,  the  one-dimensional  solutions  returned  to  their  approxi¬ 
mate  original  strain  state  after  unloading  for  interphase  yield  stresses 
below  130  MPa.  This  is  most  likely  attributed  to  the  added  matrix  plastic- 


77 


Figure  32.  Von  Mises  Stress  Contours  After  Cooldown  Using  Plane 
Stress  Solution  for  Strong  Interface  and  for  an  Interphase 
Zone  Yield  Stress  of  175  MPa 


Figure  33.  Von  Mises  Stress  Contours  After  Cooldown  Using  Plane 
y  Stress  Solution  for  Interphase  Yield  Stresses  of  125  MPa  and 

75  MPa 
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ity. 


The  effect  of  varying  the  interphase  zone  strain  hardening  parameter 
on  the  transverse  response  is  depicted  in  Figure  34.  The  one-dimensional 
solutions  again  provided  a  good  introduction  into  the  behavior  found 
using  the  finite  element  analysis.  Varying  the  interphase  zone  strain 
hardening  parameter  changes  the  secondary  slope  of  the  stress-strain 
curve  and  this  effect  can  be  easily  seen  in  the  figure. 

Unloading  is  characterized  by  elastic  followed  by  plastic  deformation, 
and  varying  the  strain  hardening  parameter  was  found  to  change  the  slope 
of  the  stress-strain  curve  during  the  plastic  part  of  unloading.  This  effect 
was  earlier  observed  in  the  one -dimensional  solutions,  but  the  ch?oge  in 
slope  was  of  opposite  sign(see  Figure  27,  page  71).  This  difference  be- 
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Figure  34.  Finite  Element  Plane  Stiess  Solutions  of  Interphase  Zone 
Strain  Hardening  Parameter  Variation 
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tween  the  one-dimensional  solutions  and  the  finite  element  results  from 
the  fact  that  the  one-dimensional  models  used  isotropic  hardening  while 
the  finite  element  models  used  kinematic  hardening.  Therefore,  plastic 
deformation  during  loading  would  increase  the  compressive  yield  stress  in 
the  one-dimensional  models  but  decrease  it  in  the  finite  element  models. 

Plane  stress  solutions  were  also  performed  on  various  interphase  zone 
sizes.  These  results  are  presented  in  Figure  35.  The  same  general  trends 
observed  in  the  one-dimensional  results  were  also  observed  here,  but  the 
effect  was  more  pronounced.  Enlarging  the  interphase  zone  lowers  the 
point  at  which  nonlinearity  first  occurs,  and  therefore  increases  the  ener¬ 
gy  dissipated  during  a  loading  and  unloading  cycle. 


Figure  35.  Finite  Element  Plane  Stress  Solutions  of  Interphase  Zone  Size 
Variation 
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Calculations  involving  the  three-dimensional  generalized  plane  strain 
model  were  costly  from  the  computer  resources  standpoint.  A  complete 
solution  would  usually  take  over  24  hours  of  processing  time,  and  many 
trial  runs  were  required  before  the  necessary  choice  of  convergence  pa¬ 
rameters  to  complete  a  solution  were  achieved.  Therefore,  only  one  set  of 
solutions  were  obtained.  This  set  consisted  of  various  interphase  zone 
yield  stresses.  These  results  are  given  in  Figure  36.  As  previously,  per¬ 
fectly  strong  and  perfectly  weak  bonded  interface  solutions  are  plotted 
alongside  for  comparison. 

The  generalized  plane  strain  solutions  demonstrated  the  same  trend  as 
the  plane  stress,  but  were  considerably  stiffer.  For  instance,  with  an  in¬ 
terphase  yield  stress  of  100  MPa  the  overall  transverse  strain  at  maximum 


Figure  36.  Finite  Element  Generalized  Plane  Strain  Solutions  of  Inter¬ 
phase  Zone  Yield  Stress  Variation  as  Compared  to  Perfectly 
Strong  and  Perfectly  Weak  Bonds 
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load  for  the  plane  stress  calculation  was  almost  double  that  found  using 
the  generalized  plane  strain  solution. 

The  additional  constraint  in  the  out-of  plane  direction  will  result  in  a 
triaxial  stress  state  throughout  the  structure.  Even  though  the  yield  stress 
of  the  interphase  zone  may  be  relatively  low,  the  triaxial  stress  state  will 
result  in  low  Von  Mises  stresses,  and  therefore  the  extent  of  yielding  will 
be  greatly  reduced. 

In  order  to  determine  the  approximate  discretization  error  associated 
with  the  finite  element  grid  chosen  for  this  study,  an  additional  grid  con¬ 
taining  approximately  four  times  the  number  of  nodes  was  developed 
(see  Figure  7,  page  37).  A  calculation  using  the  baseline  parameters  of 
Table  5  was  performed  and  a  comparison  made  with  the  conventional  grid. 
The  transverse  response,  shown  in  Figure  37,  demonstrates  excellent 
agreement  in  the  elastic  realm,  but  when  nonlinearity  begins,  the  solutions 
drift  apart.  The  greatest  discrepancy  occurs  at  maximum  load  where  the 
difference  in  overall  transverse  strain  is  0.00031,  and  remains  at  that  level 
during  the  elastic  portion  of  unloading.  The  maximum  percent  error  oc¬ 
curs  at  the  end  of  the  elastic  portion  of  unloading  and  is  12.5%.  An  error 
of  this  magnitude  may  seem  excessive,  but  it  is  important  to  remember 
that  the  solutions  represent  the  behavior  after  cooldown.  Therefore,  the 
error  is  an  accumulation  over  an  entire  thermomechanical  cycle  consisting 
of  cooldown,  transverse  loading,  and  then  unloading.  Also,  it  is  encour¬ 
aging  that  the  basic  characteristics  of  the  transverse  response  are  the 
same.  Only  a  shift  in  tne  data  has  occurred. 

In  summary,  the  finite  element  solutions  that  incorporated  an  inter¬ 
phase  zone  between  fiber  and  matrix  in  general  demonstrated  the  same 
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Figure  37.  Finite  Element  Plane  Stress  Solutions  of  Conventional  and 
Fine  Grids  Using  Baseline  Parameter  Set  of  Table  5 


trends  as  were  observed  with  the  one-dimensional  solutions.  The 
interphase  plastic  parameters(i.e.  yield  stress  and  strain  hardening)  were 
found  to  have  the  greatest  effect  on  the  overall  behavior.  Increasing  the 
yield  stress  resulted  in  increasing  the  point  of  nonlinearity,  and  increasing 
the  strain  hardening  resulted  in  increasing  the  secondary  slope.  Most  im¬ 
portantly,  modeling  the  interface  as  a  third  phase(interphase)  elastic-plas¬ 
tic  zone  allowed  for  continuous  variation  in  transverse  behavior  from 
strong  to  weak  bonding(see  Figure  31  on  page  76).  Thus,  such  a  model 


allows  for  additional  constraints  at  the  interface  region  that  may  be 


tailored  to  match  a  particular  composite’s  characteristics.  Other 
attributes  of  the  interphase  zone  model  included  significant  energy  dissi- 
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pation  during  a  loading  and  unloading  sequence,  and  a  relaxation  of  the 
thermal  residual  stresses  due  to  the  interphase  zone  yielding  during 
cooldown. 

The  main  difference  between  the  finite  element  and  one-dimensional 
solutions  was  that  the  finite  element  solutions  always  experienced  perma¬ 
nent  strain  after  a  loading  and  unloading  cycle.  This  was  most  likely  due 
to  the  added  plasticity  of  the  matrix  in  the  finite  element  analysis. 

Figure  38  gives  a  transveru-  response  comparison  between  the  one-di¬ 
mensional  and  plane  stress  finite  element  analysis  solutions  of  the  base¬ 
line  parameter  set.  The  one-dimensional  solution  does  an  excellent  job  of 
predicting  the  general  characteristics,  pi  iducing  only  a  12%  error  in  total 
strain  at  maximum  load. 
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Figure  38.  Comparison  of  One-Dimensional  and  Finite  Element 

Transverse  Response  for  Baseline  Parameter  Set  of  Table  4 
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A  comparison  of  the  nonlinear  characteristics  as  calculated  by  the  one¬ 
dimensional  and  finite  element  analyses  is  depicted  in  Figure  39.  The 
one-dimensional  analysis  tended  to  overpredict  the  point  of  nonlinearity 
and  underpredict  the  secondary  slope.  The  point  of  nonlinearity  varied 
with  interphase  yield  stress  in  a  bilinear  nature.  This  behavior  was  more 
pronounced  in  the  one-dimensional  solutions,  and  it  is  interesting  to  note 
that  the  bilinearity  begins  at  an  interphase  yield  stress  of  approximately 
125  MPa.  This  is  the  same  yield  stress  at  which  a  relaxation  in  the  ther¬ 
mal  residual  stresses  during  cooldown  began  to  occur  in  the  one-dimen¬ 
sional  solutions.  However,  as  a  final  note,  considering  the  simplifying  as¬ 
sumptions  made  in  the  one-dimensional  analysis,  the  overall  agreement 
between  it  and  the  finite  element  calculations  is  excellent. 

4.3  Results  of  Combined  Interphase  Zone  /  Weak  Bond  Models 

The  previous  section  points  out  that  a  major  difficulty  in  modeling  the 
fiber/matrix  interface  as  a  third  zone  possessing  independent  properties  is 
the  relaxation  of  the  thermal  residual  stresses  during  cooldown  as  the  in¬ 
terphase  zone  yields.  Another  drawback  is  that  all  models  considered  so 
far  have  assumed  that  the  conditions  are  identical  at  all  points  along  the 
interface  (i.e.  the  interface  behaves  the  same  everywhere).  Therefore,  this 
section  will  examine  models  possessing  an  interface  with  combined  ef¬ 
fects.  Part  of  the  interface  will  be  modeled  as  perfectly  weak  and  part  as 
possessing  an  interphase  zone.  The  development  of  these  models  is  dis¬ 
cussed  in  section  3.5.  The  energy  dissipation  capability  of  an  interphase 
zone  coupled  with  a  perfectly  weak  bond  constraint  to  prevent  the  relax¬ 
ation  of  thermal  stresses  could  provide  a  superior  micromechanics  model. 
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Figure  39.  Comparison  of  Nonlinear  Characteristics  Between  the  One- 
Dimensional  and  Finite  Element  Solutions 
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4.3.1  One-Dimensional  Results  of  a  Combined  Interphase  Zone  / 
Weak  Bond  Model 


The  loading  sequence  will  be  the  same  as  in  section  4.2,  cooldown  to 
room  temperature  followed  by  transverse  normal  loading  and  unloading. 
The  baseline  material  properties  and  model  dimensions  will  be  the  same 
as  listed  in  Table  4.  An  additional  dimension  not  listed  in  Table  4  is  re¬ 
quired.  This  dimension,  d,  shown  in  Figure  9  is  a  representation  of  the 
portion  of  the  interface  modeled  as  perfectly  weak  and  will  be  given  ref¬ 
erenced  to  the  fiber  size,  a.  Therefore,  since  shear  effects  are  neglected, 
the  ratio  d/a  provides  an  approximate  representation  of  the  percentage  of 
the  interface  modeled  as  perfectly  weak  and  will  be  referred  to  as  the 
weak  bond  ratio.  Since  it  was  found  in  section  4.2  that  the  plastic  charac¬ 
teristics  of  the  interphase  zone  effected  composite  behavior  much  more 
than  the  other  interphase  zone  parameters,  the  parameters  varied  in  this 
portion  of  the  study  were  only  interphase  zone  yield  stress,  strain  harden¬ 
ing,  and  the  weak  bond  ratio.  For  validation  of  the  combined  interphase 
zone  /  weak  bond  one-dimensional  solutions  the  weak  bond  ratio  was  set 
to  both  0%  and  100%  to  see  if  the  behavior  would  approach  that  of  an  all 
interphase  zone  model  and  a  perfectly  weak  bond  interface  model. 

Figure  40  depicts  the  transverse  loading  and  unloading  behavior  for  these 
two  cases  as  well  as  a  weak  bond  ratio  of  50%.  A  weak  bond  ratio  of  zero 
matches  exactly  the  results  for  the  one-dimensional  model  possessing  only 
an  interphase  zone.  However,  there  is  some  energy  dissipation  that  occurs 


for  a  weak  bond  ratio  of  100%  that  is  not  seen  in  a  perfectly  weak  bond 


model.  This  discrepancy  results  from  the  fact  that  in  the  one-dimensional 
solutions  a  portion  of  the  interphase  zone  above  the  fiber  remains  intact  in 
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Figure  40.  One-Dimensional  Solutions  of  Combined  Interphase  Zone  / 
Weak  Bond  for  Weak  Bond  Ratios  from  0.0  to  1.0 


spite  of  a  100%  weak  bond  ratio.  The  plastic  deformation  of  this  remain¬ 
ing  third  phase  material  causes  the  energy  dissipation. 

A  rather  surprising  result  from  Figure  40  is  the  behavior  of  the  50% 
weak  bond  ratio.  The  behavior  of  this  curve  is  difficult  to  interpret  with¬ 
out  examining  the  stresses  and  displacements  in  each  region.  A  summary 
of  what  occurs  during  loading  and  unloading  is  given  in  Figure  41.  The 
curve  may  be  divided  into  several  distinct  linear  fields.  These  character¬ 


istic  linear  fields  existed  for  all  the  one-dimensional  combined  interphase 


zone  /  weak  bond  solutions.  Only  the  slope  and  size  of  the  fields 
changed. 
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Figure  41.  One-Dimensional  Solution  of  Combined  Interphase  Zone  / 
Weak  Bond  With  Loading  Sequence  Explanations 


An  investigation  into  the  effect  of  the  weak  bond  ratio  on  the  trans¬ 
verse  behavior  is  depicted  in  Figure  42.  A  few  characteristics  of  the  one¬ 
dimensional  model  when  increasing  the  weak  bond  ratio  are  worth  men¬ 
tioning.  First,  the  stress  at  which  fiber/matrix  separation  occurs  increas¬ 
es.  Also,  the  slope  in  the  second  linear  field(after  separation  but  before 
interphase  zone  yielding)  decreases  and  approaches  the  slope  that  exists 
after  interphase  zone  yield.  The  combined  model  was  found  to  prevent  re 
laxation  of  the  thermal  residual  stresses  by  keeping  the  interphase  zone 
from  yielding  during  cooldown.  Therefore,  as  the  weak  bond  ratio  was  in 
creased  the  solutions  approached  that  of  a  completely  weak  bonded  inter¬ 
face  but  still  retained  an  energy  dissipation  mechanism  in  the  form  of  a 
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Figure  42.  One-Dimensional  Solutions  of  Combined  Interphase  Zone  / 
Weak  Bond  for  Various  Weak  Bond  Ratios 


small  portion  of  the  interface  modeled  as  a  third  phase  material. 

Since  the  interphase  zone  yield  stress  had  a  pronounced  effect  on 
the  one-dimensional  model  of  section  4.2.1,  its  effect  on  the  combined 


model  was  also  examined.  Figure  43  presents  an  interphase  yield  stress 
variation  for  a  weak  bond  ratio  of  50%.  As  expected,  an  increase  in  yield 
stress  lengthens  the  portion  of  the  curve  between  fiber/matrix  separation 
and  interphase  yielding  but  shortens  the  one  between  interphase  yield  and 
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stress  variation  at  a  higher  weak  bond  ratio  as  depicted  in  Figure  44.  The 
same  effect  is  occurring  but  it  is  manifesting  itself  in  what  seems  to  be  a 
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Figure  43.  One-Dimensional  Solutions  of  Combined  Interphase  Zone  / 
Weak  Bond  With  a  Weak  Bond  Ratio  of  0.50  for  Various  In¬ 
terphase  Zone  Yield  Stresses 


change  in  slope  of  the  upper  portion  of  the  curve. 

The  final  investigation  where  a  one-dimensional  solution  was 
employed  involved  the  effect  of  interphase  zone  strain  hardening  for  the 
combined  model.  This  is  shown  in  F  gure  45,  and  as  expected  an  increase 
in  strain  hardening  results  in  increasing  slope  during  the  portion  of  load¬ 
ing  where  the  interphase  zone  yields. 

For  a  one-dimensional  model  the  portions  of  the  loading  curve  that  are 
controlled  by  the  various  mechanisms  described  in  Figure  41  are  well  de¬ 
fined,  but  for  a  two  or  three  dimensional  analysis  where  the  fiber  is  circu¬ 
lar  rather  than  square  this  is  not  the  case.  Not  all  of  the  interphase  zone 
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Figure  44.  One-Dimensional  Solutions  of  Combined  Interphase  Zone  / 
Weak  Bond  With  a  Weak  Bond  Ratio  of  0.90  for  Various  In¬ 
terphase  Zone  Yield  Stresses 


yields  at  the  same  transverse  load.  After  initial  yield  of  a  small  area  of 
the  interphase  any  further  load  will  cause  a  larger  area  to  yield  and  so  on 
until  the  entire  region  is  plastic.  In  addition,  a  simplified  one-dimension¬ 
al  model  may  provide  a  general  understanding  of  how  varying  a  specific 
parameter  will  induce  trends  in  the  data,  but  the  calculations  are  far  from 
an  exact  solution.  Therefore,  an  investigation  into  a  combined  interphase 
zone  /  weak  bond  micromechanics  model  using  more  exact  methods  must 
be  sought. 
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Figure  45.  One-Dimensional  Solutions  of  Combined  Interphase  Zone  / 
Weak  Bond  for  Various  Strain  Hardening  Parameters 


4.3.2  Finite  Element  Results  of  a  Combined  Interphase  Zone  /  Weak 
Bond  Model 

A  description  of  the  finite  element  models  used  to  analyze  a  combined 
interphase  zone  /  weak  bond  interface  between  fiber  and  matrix  is  given  in 
section  3.5.  As  mentioned  in  the  development  of  the  finite  element  mod¬ 
els,  a  choice  as  to  the  positioning  of  the  interphase  zone  along  the 
fiber/matrix  interface  was  necessary.  A  position  was  chosen  such  that  the 
interphase  zone  was  centered  symmetrically  about  a  45°diagonal  of  the 
representative  volume  element(see  Figure  10,  page  43). 

Furthermore,  the  weak  bond  ratio  for  the  finite  element  calculations 
was  defined  to  be  the  length  along  the  interface  modeled  as  perfectly 
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weak  divided  by  the  total  interface  length.  The  properties  used  were 
those  listed  in  Table  5  on  page  75,  and  both  the  weak  bond  ratio  and  inter¬ 
phase  zone  yield  stress  were  varied  to  determine  their  effect  on  the  over¬ 
all  composite  characteristics.  Also,  the  cooldown  and  loading  sequence 
was  identical  to  that  used  for  the  one-dimensional  solutions. 

Since  several  weak  bond  ratios  were  analyzed,  the  solutions  are  pre¬ 
sented  in  separate  figures  to  avoid  excessive  clutter.  The  results  for  weak 
bond  ratios  of  50%  and  higher  are  presented  in  Figure  46,  and  the  results 
for  weak  bond  ratios  lower  than  50%  are  given  in  Figure  47.  One  of  the 
most  noticeable  characteristics  is  the  tremendous  difference  in  transverse 
response  between  the  perfectly  weak  (100%  weak  bond  ratio)  and  the  87% 
weak  bond  ratio.  Even  when  a  relatively  small  section  between  the  fiber 
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Figure  46.  Finite  Element  Plane  Stress  Solutions  of  Combined  Inter¬ 
phase  Zone  /  Weak  Bond  Models  for  Various  Weak  Bond 
Ratios  -  High  Range 
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Figure  47.  Finite  Element  Plane  Stress  Solutions  of  Combined  Inter¬ 
phase  Zone  /  Weak  Bond  Models  for  Various  Weak  Bond  Ra¬ 
tios  -  Low  Range 


and  matrix  is  modeled  as  an  interphase  zone,  the  structure  becomes  much 
stiffer.  The  one-dimensional  solutions  presented  a  gradual  change  to  a 
perfectly  weak  bond  as  the  weak  bond  ratio  was  varied.  This  discrepancy 
in  the  one-dimensional  solutions  is  possibly  due  to  neglecting  constraints 
on  the  interphase  zone  from  the  surrounding  material.  In  the  finite  ele¬ 
ment  model  the  interphase  zone  is  entirely  encased  in  either  matrix  or 
fiber  mateiial.  Therefore,  any  deformation  of  this  zone  will  necessarily 
deform  the  surrounding  material,  but  the  fiber  and  matrix  are  considerably 
stronger  and  will  resist  it.  Also,  in  the  finite  element  solutions  varying 
the  weak  bond  ratio  does  not  effect  the  point  of  fiber/matrix  separation  or 
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the  permanent  strain  existent  after  unloading.  However,  the  general  trend 
seen  in  the  upper  portion  of  the  transverse  response  curve  is  very  similar 
to  the  one-dimensional  solutions  since  at  the  higher  load  levels  a  weaken¬ 
ing  of  the  structure  occurs  as  the  weak  bond  ratio  is  increased. 

The  low  range  of  weak  bond  ratios  (below  50%)  also  demonstrated 
surprising  behavior  when  the  weak  bond  ratio  was  changed  from  0%  to 
7%.  Again  a  jump  in  the  results  occurred  as  the  model  was  changed  from 
taking  the  entire  interface  as  an  interphase  zone  to  a  combined  interphase 
zone  /  weak  bond  model.  The  same  possible  explanation  applies  here  as 
was  given  for  the  87%  to  100%  weak  bond  ratio  jump.  The  combined  in¬ 
terphase  zone  /  weak  bond  models  require  the  interphase  to  be  completely 
surrounded  by  fiber  and  matrix  material,  creating  a  constraint  on  the  in¬ 
terphase  boundary.  Hence,  an  interesting  development  has  occurred.  The 
combined  model  is  stiffer  than  either  of  the  two  by  themselves.  Also,  for 
all  weak  bond  ratios  greater  than  0%  relaxation  of  the  thermal  residual 
stresses  during  cooldown  did  not  occur.  As  a  final  note  on  varying  the 
weak  bond  ratio,  very  little  effect  is  seen  on  the  transverse  response  for 
weak  bond  ratios  below  50%. 

An  examination  into  how  varying  the  interphase  zone  yield  stress  ef¬ 
fects  the  combined  model  is  given  in  Figure  48  for  a  weak  bond  ratio  of 
87%.  The  results  demonstrate  how  increasing  the  interphase  yield  stress 
stiffens  the  structure  and  decreases  the  energy  dissipated  per  cycle.  This 
same  effect  was  observed  in  the  one-dimensional  solutions.  As  the  yield 
stress  was  increased  in  the  finite  element  solutions  the  curve  approached 
nonlinear  elastic  behavior  where  the  nonlinearity  results  from  separation 
of  the  portion  of  the  interface  modeled  as  a  weak  bond. 
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Figure  48.  Finite  Element  Plane  Stress  Solutions  of  Combined  Inter¬ 
phase  Zone  /  Weak  Bond  With  a  Weak  Bond  Ratio  of  87%  for 
Various  Interphase  Zone  Yield  Stresses 


Two  solutions  using  the  three-dimensional  generalized  plane  strain 
model  with  a  combined  interphase  zone  /  weak  bond  were  obtained.  They 
were  at  weak  bond  ratios  of  7%  and  87%.  Figure  49  presents  the  7%  re¬ 
sults  as  compared  to  the  plane  stress  solution.  Figure  50  presents  the  87% 
results  in  the  same  manner.  As  in  previous  generalized  plane  strain  solu¬ 
tions  from  earlier  sections,  the  calculations  result  in  a  stiff er  overall 
structure  and  decreased  energy  dissipation  per  cycle.  However,  in  the  lin¬ 
ear  elastic  region  before  fiber/matrix  separation  the  two  types  of  analysis 
agree  very  well.  Also,  in  the  87%  weak  bond  ratio  solution  an  increase  in 
the  point  of  fiber/matrix  separation  of  approximately  50  MPa  occurs  for 
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Figure  49.  Finite  Element  Generalized  Plane  Strain  Solution  of 

Combined  Interphase  Zone  /  Weak  Bond  at  7%  Weak  Bond 
Ratio  as  Compared  to  Plane  Stress  Solution 

the  generalized  plane  strain  analysis  over  the  plane  stress. 

In  summary,  the  combined  interphase  zone  /  weak  bond  models  have 
demonstrated  the  ability  to  model  energy  dissipation  at  the  interface,  and 
at  the  same  time  prevent  the  relaxation  of  thermal  residual  stresses  which 
was  a  drawback  with  the  models  that  used  only  an  interphase  zone.  How¬ 
ever,  the  combined  model  offers  little  flexibility  in  the  transverse  behav¬ 
ior  because  it  greatly  stiffens  the  solutions  over  both  the  perfectly  weak 
interface  and  all  interphase  zone  models.  Also,  for  all  finite  the  element 
solutions,  permanent  strain  still  exists  after  unloading  which  has  not  been 
observed  in  experiment.  Finally,  the  one-dimensional  analysis  deviated 
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Figure  50.  Finite  Element  Generalized  Plane  Strain  Solution  of 

Combined  Interphase  Zone  /  Weak  Bond  at  87%  Weak  Bond 
Ratio  as  Compared  to  Plane  Stress  Solution 


more  from  the  finite  element  solutions  in  this  section  than  was  observed 
in  section  4.2.  The  simplifying  assumptions  of  the  one-dimensional  anal¬ 
ysis  were  not  as  capable  of  handling  the  more  complex  combined  inter¬ 
phase  zone  /  weak  bond  model  of  the  interface. 
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V.  Conclusions 


The  response  of  a  micromechanical  model  is  highly  dependent  on  the 
fiber/matrix  interface  condition.  An  analyst  must  understand  the  proper 
interface  conditions  required  for  a  particular  composite  before  undertak¬ 
ing  an  investigation.  This  study  has  examined  various  interface  condi¬ 
tions  to  determine  their  effect  on  the  overall  transverse  normal  character¬ 
istics  using  both  a  simplified  one-dimensional  analysis  and  the  finite  ele¬ 
ment  method. 

Traditionally,  a  perfectly  strong  interfacial  bond  is  assumed  for  poly¬ 
mer  or  epoxy  based  composites.  However,  this  assumption  is  erroneous 
for  high  temperature  titanium  based  metal  matrix  composites.  Experi¬ 
ments  indicate  that  the  interfacial  bond  for  these  composites  fails  early 
during  initial  transverse  loading [33].  In  addition,  the  metal  matrix  typical¬ 
ly  possesses  a  much  higher  thermal  coefficient  of  expansion  than  the 
fiber.  Therefore,  relatively  large  thermal  residual  stresses  are  present  in 
the  composite  after  cooldown  from  the  processing  temperature,  and  final¬ 
ly,  a  reaction  zone  is  present  at  the  fiber/matrix  interface  where  the  fiber 
has  penetrated  into  and  reacted  with  the  matrix  forming  a  third  phase  ma¬ 
terial  which  is  likely  to  possess  its  own  unique  properties. 

Large  thermal  residual  compressive  stresses  exist  at  the  interface  after 
cooldown  in  a  high  temperature  metal  matrix  composite.  These  compres¬ 
sive  stresses  must  be  overcome  before  the  matrix  will  separate  from  the 
fiber.  Previous  researchers  have  successfully  modeled  this  phenomena  by 
using  a  perfectly  weak  interface^;  3].  However,  difficulties  arise  from 
such  models  due  to  their  lack  of  sufficient  constraint  at  the  interface. 
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Their  overall  transverse  response  cannot  be  modified  to  model  a  partial 
bonding  at  the  interface.  As  a  result,  the  weak  bond  solutions  do  not  al¬ 
ways  agree  with  experiment(see  Figure  2,  page  16). 

Hence,  parametric  studies  into  how  various  interface  conditions  effect 
the  overall  composite  characteristics  will  provide  insight  into  promising 
methods  of  modeling  the  complex  fiber/matrix  interface  region.  This 
study  has  examined  both  strong  and  weak  bond  models,  investigating 
fiber/matrix  separation  and  slip.  In  addition,  an  interphase  zone  between 
fiber  and  matrix  has  been  added  in  an  effort  to  model  the  energy  dissipa¬ 
tion  mechanisms  and  additional  constraints  that  cannot  be  modeled  using 
a  weak  bond.  Finally,  a  combined  interphase  zone  /  weak  bond  interface 
model  was  examined. 

An  important  conclusion  from  the  weak  bond  model  is  that  not  only 
fiber/matrix  separation  but  also  slip  along  the  interface  was  found  to  have 
a  pronounced  effect  "n  the  overall  characteristics.  If  friction  is  absent  at 
the  interface,  then  the  initial  modulus  that  occurs  before  fiber/matrix  sep¬ 
aration  is  considerably  less  than  for  a  strong  bond.  This  is  due  to  interfa¬ 
cial  slippage.  When  friction  is  added,  slip  is  arrested  and  the  initial  mod¬ 
ulus  increases  and  approaches  that  of  a  strong  bond.  However,  after 
fiber/matrix  separation  the  added  constraint  from  friction  vanishes  since 
the  fiber  and  matrix  are  no  longer  in  contact,  and  hence,  the  secondary 
slope  after  separation  remains  unchanged  from  an  interface  without  fric- 
tion(see  Figure  17,  page  54). 

If  the  fracture  surface  between  fiber  and  matrix  is  jagged  in  nature, 
then  separation  and  slip  will  be  independent.  In  other  words,  if  the  sepa¬ 
ration  distance  between  fracture  surfaces  is  small,  slip  will  still  not  be  al- 
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lowed.  Adding  such  a  constraint  that  prevents  slip  throughout  the  load 
range  was  found  to  greatly  stiffen  the  structure.  In  fact,  the  added  con¬ 
straint  stiffened  the  overall  structure  so  much  that  the  stresses  in  the  ma¬ 
trix  remained  elastic,  and  therefore,  a  nonlinear  elastic  curve  resulted. 

The  transverse  response  with  and  without  slip  curves  were  shown  to 
bound  previously  published  experimental  data(see  Figure  24,  page  64). 

Modeling  the  nonlinearities  at  the  interface  by  assuming  an  interphase 
zone  possessing  unique  material  properties  allowed  for  excellent  flexibili¬ 
ty  for  including  additional  constraints  not  possible  with  the  weak  bond 
model.  Also,  the  motivation  for  incorporating  an  interphase  zone  model 
is  physically  realistic  in  that  it  is  based  on  experimental  observations  of  a 
reaction  zone  between  the  fiber  and  matrix[39].  Cracks  and  other  imper¬ 
fections  present  in  the  reaction  zone  are  modeled  through  the  elastic-plas¬ 
tic  properties  of  the  interphase.  The  plastic  properties  (i.e.  yield  stress 
and  strain  hardening)  of  the  interphase  zone  were  found  to  have  the  great¬ 
est  effect  on  the  transverse  characteristics  of  the  composite.  The  zone 
size  also  had  considerable  impact  on  the  results,  but  the  elastic  properties 
of  the  interphase  such  as  Young’s  modulus  and  thermal  coefficient  of  ex¬ 
pansion  exhibited  only  slight  influence.  In  addition  to  possessing  the  ca¬ 
pability  for  modeling  a  partially  bonded  interface,  an  advantage  in  using 
an  interphase  zone  was  that  it  allowed  for  significant  energy  dissipation 
during  a  loading  and  unloading  sequence.  Finally,  the  interphase  zone 
models  were  found  to  produce  a  relaxation  of  the  thermal  residual  stresses 
as  the  interphase  exceeded  its  yield  stress  in  compression  during 
cooldown(see  Figure  33  on  page  79).  However,  for  most  materials  the 
initial  yield  stress  in  compression  is  greater  than  in  tension,  but  the 
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analysis  methods  used  in  this  study  did  not  have  the  capability  to  specify 
initial  yield  stresses  that  were  different  in  compression  and  tension. 

The  combined  interphase  zone  /  weak  bond  interface  was  used  for 
modeling  a  portion  of  the  interface  as  if  a  perfectly  weak  bond  existed  and 
the  remainder  as  if  an  interphase  zone  existed  between  fiber  and  matrix. 

It  offered  an  advantage  over  the  all  interphase  zone  model  in  that  it  al¬ 
lowed  for  energy  dissipation  during  a  loading  and  unloading  sequence 
while  at  the  same  time  preventing  a  relaxation  of  the  thermal  residual 
stresses.  A  drawback  of  the  combined  interphase  zone  /  weak  bond  inter¬ 
face  was  that  for  the  finite  element  analysis  it  did  not  allow  for  flexibility 
in  the  transverse  response  characteristics.  For  instance,  the  combined 
model  was  in  general  substantially  stiffer  than  either  method  used  alone, 
and  varying  the  parameters,  such  as  the  interphase  material  properties  and 
the  weak  bond  ratio,  had  little  impact  on  the  transverse  response.  This 
behavior  is  most  likely  due  to  the  fact  that  the  interphase  is  completely 
enclosed  in  either  fiber  or  matrix,  and  hence,  any  deformation  of  its 
boundary  will  require  equal  deformation  of  the  surrounding  material 
which  is  considerably  stronger.  However,  the  one-dimensional  solutions 
did  not  experience  this  difficulty  since  the  simplifying  assumptions  re¬ 
moved  all  but  a  single  displacement  constraint.  A  summary  of  the  various 
methods  for  modeling  the  interface  and  their  characteristics  are  presented 
in  Table  6. 

In  review  of  the  analysis  methods  used  in  this  study,  both  the  finite 
element  and  simplified  one-dimensional  results  are  presented  with  a  high 
degree  of  confidence.  This  confidence  comes  from  agreement  between  the 
methods  themselves  and  validation  from  a  previous  researcher’s  calcula- 
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Table  6.  Summary  of  Interface  Models  and  Their  Characteristics 


Interface  Model  Characteristics 


Strong 

-  Linear  elastic  to  matrix  yield 

-  Sufficient  for  determining  composite  behavior 
before  fiber/matrix  separation 

Weak 

-  Captures  fiber/matrix  separation. 

-  Energy  dissipation  in  the  form  of  matrix  plastic  flow. 
Results  in  permanent  strain  after  unloading. 

-  Change  in  interfacial  friction  has  little  effect  on 
secondary  slope  of  transverse  response  curve 

Weak  with 

fiber/matrix 

slip 

constrained 

-  Captures  fiber/matrix  separation. 

-  Greatly  stiffens  transverse  response  so  that  no  matrix 
plastic  flow  occurs  even  under  substantial  loading. 

-  Results  in  nonlinear  elastic  response. 

Interphase 

Zone 

-  Simulates  fiber/matrix  separation  in  the  form  of 
interphase  yielding. 

-  Capable  of  modeling  partial  bonding  by  varying  the  in¬ 

terphase  zone  plastic  properties. 

-  Allows  for  energy  dissipation  in  the  form  of  both 
interphase  and  matix  plastic  flow. 

-  Relaxation  of  thermal  residual  stresses  occurs  during 
cooldown. 

-  Permanent  strain  exists  after  unloading,  but  not  as 
significant  as  with  the  perfectly  weak  bond  model. 

Combined 
Interphase 
Zone  /  Weak 
Bond 

-  Finite  element  solutions  produce  stiffer  transverse 
response  than  observed  using  either  interphase  zone 
or  weak  bond  models  alone,  and  parameter  variations 
have  little  effect  on  it. 

-  Energy  dissipation  in  the  form  of  both  interphase 
and  matrix  plastic  flow. 

-  Prevents  relaxation  of  thermal  stresses  that  was  expe¬ 
rienced  in  the  interphase  zone  model. 

-  Permanent  strain  after  unloading  similar  to  interphase 
zone  model. 
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tions(see  Figure  18,  page  55). 

A  surprising  facet  of  the  solution  methods  was  that  a  plane  stress  rath¬ 
er  than  a  plane  strain  analysis  was  found  to  more  closely  approximate  a 
generalized  plane  strain  analysis.  This  was  concluded  during  the  linear 
elastic  portion  of  the  study  when  using  a  strong  interfacial  bond.  The  ex¬ 
planation  for  this  apparent  anomaly  is  thought  to  be  linked  to  the  large 
thermal  strains  which  are  an  order  of  magnitude  greater  than  the  mechani¬ 
cal  strains.  The  plane  strain  analysis  prevents  any  thermal  contraction  in 
the  out-of-plane  direction.  On  the  other  hand,  both  th  ;  plane  stress  and 
generalized  plane  strain  analyses  will  allow  out-of -plane  thermal  contrac¬ 
tion. 

The  one-dimensional  solutions  were  in  general  very  good  considering 
the  simplifying  assumptions.  The  one-dimensional  results  from  the  inter¬ 
phase  zone  model  demonstrated  excellent  agreement  with  the  finite  ele¬ 
ment  solutions.  All  trends  due  to  a  particular  parameter  variation  were  in 
accord,  and  for  the  most  part,  even  the  numerical  values  were  close  to  the 
same.  In  fact,  the  level  of  sophistication  in  the  one-dimensional  model 
was  greater  than  required.  The  interphase  region  above  the  fiber  was 
found  to  have  only  a  slight  effect  on  the  solutions,  and  could  have  been 
removed  without  greatly  effecting  the  results. 

The  finite  element  code  MSC/NASTRAN  was  exceptional  in  solving 
the  nonlinear  material  problem.  As  mentioned  earlier,  the  solutions  dem¬ 
onstrated  excellent  agreement  with  previous  calculations  presented  in  the 
literature^;  3].  Also,  a  converged  solution  at  all  load  points  was  generally 
achieved  after  no  more  than  three  NASTRAN  job  submittals.  Hence,  the 
analysis  methods  used  in  this  study  were  found  to  provide  excellent  in- 
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sight  into  micromechanical  analysis  problems  that  possess  nonlinear  ter- 
face  conditions.  As  demonstrated  'n  this  study,  the  solution  of  such  prob¬ 
lems  are  essential  when  dealing  with  imperfect  interfacial  bonds  in  the 
presence  of  thermal  residual  stresses. 
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VI.  Recommendations  and  Suggestions 


An  interphase  zone  between  the  fiber  and  matrix  was  incorprated  in  an 
attempt  to  model  a  fiber/matrix  bond  other  than  perfectly  strong  or  weak, 
in  effect,  a  partial  bond.  Rather  than  modeling  a  partial  bond  in  such  a 
fashion,  one  might  consider  developing  a  new  type  of  element  similar  to 
the  gap  element.  However,  it  would  need  to  possess  additional  constraint 
capabilities,  such  as  transverse  and  longitudinal  stiffness  existing  both 
before  and  after  separation  and  a  friction  capability  in  the  longitudinal 
direction.  Developing  such  an  element  might  be  a  formidable  task,  and 
the  additional  capabilities  would  have  limited  use  for  other  types  of 
applications,  but  nevertheless,  such  an  element  could  be  tailored  to  model 
any  level  of  interfacial  bonding  from  perfectly  strong  to  perfectly  weak. 

Also,  since  the  types  of  composites  considered  in  this  analysis  are  en¬ 
visioned  for  high  temperature  applications,  studies  of  the  interfacial  ef¬ 
fects  during  thermomechanical  loading  cycles  where  both  the  temperature 
and  applied  load  are  changing  simultaneously  would  be  highly  beneficial. 
The  results  presented  in  this  paper  are  from  applied  transverse  loading  at 
constant  temperature.  An  aircraft  engine  component  such  as  a  turbine 
blade  will  undergo  extensive  load  and  temperature  variations  from  engine 
start  up  to  shut  down.  The  composite  stress-strain  history  during  such  a 
thermomechanical  cycle  is  instrumental  in  determining  the  useful  life  of 
the  component.  Therefore,  it  is  imperative  to  understand  how  the  integri- 
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thermomechanical  loading  cycles. 
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Finally,  a  few  suggestions  are  in  order  concerning  the  use  of 
MSC/NASTRAN  for  solving  the  nonlinear  problems  presented  in  this 
study.  For  solutions  involving  the  use  01  gap  elements  it  is  recommended 
to  update  the  stiffness  matrices  after  each  iteration.  The  author  found  that 
incorrect  results  would  occur  if  this  recommendation  was  not  followed. 
Also,  when  using  the  elastic-plastic  interphase  zone  micromechanics 
model,  it  is  recommended  to  set  the  values  of  Young’s  modulus  and  ther¬ 
mal  coefficient  of  expansion  of  the  interphase  zone  as  close  to  those  of 
the  matrix  properties  as  is  practical.  Otherwise,  convergence  may  be  dif¬ 
ficult  to  achieve. 
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Appendix  A:  Sample  MSC/NASTRAN  Data  Deck 


ID  DROBERTS , TNOINTER 

TIME  120 

SOL  66 

DIAG  50 

ALTER  73  $ 

OUTPUT2  CSTM, GPL, GPDT, EPT, MPT//0/11 
OUTPUT2  GEOM2 , GEOM3 , GEOM4 , , //0/11  $ 
ALTER  1036  $ 

CUTPUT2  OUGV1, OES1, OSTR1, , //0/11  $ 
OUTPUT2  ,  ,  ,  , // -9/11  $ 

CEMD 

TITLE  =  NO  INTER,  WEAK,  PL  STRESS 
DISPLACEMENT  =  ALL 
STRESS  =  ALL 
SHALL  =  ALL 
SUPER  =  ALL 

SPC  =  1 

MPC  =  5 

SUBCASE  1 

SUBTITLE  =  TEMP  970 
TEMP (LOAD)  =  1 
NLP ARM  =  1 

SUBCASE  2 

TEMr (LOAD;  =  10 
NLPARM  =  2 

SUBCASE  3 

SUBTITLE  =  C100+TEMP  ,  100+TEMP 
TEMF (LOAD)  =  10 
LOAD  =  2 

NLPARM  =  3 

SUBCASE  4 

SUBTITLE  =  C150+TEMP  ,  150+TEMP 
TEMP (LOAD)  =  10 
LOAD  =  3 

NLPARM  =  3 

SUBCASE  5 

SUBTITLE  =  C175+TEMP  ,  175+TEMP 
TEMP (LOAD)  =  10 
LOAD  =  4 

NLPARM  =  3 

SUBCASE  6 

SUBTITLE  =  C200+TEMP  ,  200+TEMF 
TEMP (LOAD)  =  10 
LOAD  =  5 

NLPARM  =  3 

SUBCASE  7 

SUBTITLE  =  C225+TEMP  ,  225+TEMP 
TEMP (LOAD)  =  10 


Executive  Control 
Deck 


Case  Control  Deck 
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LOAD  =  6 

NLPARM  =  3 

SUBCASE  8 

SUBTITLE  =  C250+TEMP  ,  250+TEMP 
TEMP  (LOAD)  =  10 
LOAD  =  7 

NLPARM  =  3 

BEGIN  BULK 


NLPARM  1 

1 

ITER 

1 

-40  PW 

NO 

NLPARM  2 

8 

ITER 

1 

-20  PW 

NO 

NLPARM  3 

5 

ITER 

1 

-20  PW 

NO 

GRID* [16  spaces]  291 
*GA  1 

0. 

0 

.0148493 

0 

. 0148493+GA 

1 

GRID*[16 spaces]  292 
*GA  2 

0. 

0 

.0128879 

0 

.  0165803+GA 

2 

GRID* [16 spaces]  293 
*GA  3 

0. 

0 

.0108472 

0 

. 017981702+GA 

3 

[grid  points  continued] 


CQUAD4 

545 

1 

291 

300 

301 

292 

0. 

CQUAD4 

546 

1 

292 

301 

302 

293 

0. 

CQUAD4 

547 

1 

293 

302 

303 

294 

0. 

CQUAD4 

548 

1 

294 

303 

304 

295 

0. 

CQUAD4 

549 

1 

295 

304 

305 

296 

0. 

[quad4 

i  elements  continued] 

I 


CTRIA3 

686 

1  296  297  444 

0 

CTP.IA3 

687 

1  297  298  444 

[tria3  elements  continued] 

0 

CGAP 

817 

2 

453 

371 

1 

CGAP 

818 

2 

452 

370 

2 

CGAP 

819 

2 

451 

369 

3 

[gap  elements  continued] 


CORD2R*[16  spaces]  1 

0 

6 . 208301E- 08 

*COIl  [17  spaces]  0. 

6 . 208301E- 08 

0.0699999 

*C012 

0. 

0. 

CORD2R*[  16  spaces]  2 

0 

0.00490102 

*C021  [17  spaces]  0. 

0.00490102 

0.0698274 

*0022 

0. 

0. 

0 , 0699999+CO11 
0 . 01+CO12 

0. 

0.0698274+CO21 
0 . 01+CO22 


0. 


[coordinate  systems  definitions  continued] 

I 

I 
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MATlJll spaces]!  400000.153846.2 

. 37 . 82E- 09  5.E-06  1000. 

+MA 

1 

+MA  1 

1500.  1500. 

68. 

MATlfll  spaces]2  100000.38461. 

54 

. 37 . 82E-09  .00001  1000. 

+MA 

2 

+MA  2 

1500.  1500. 

68. 

PSHELL* 

[16  spaces]  1 

1  .01 

1+PA 

1 

*PA  1 

1. 

1  .8333333 

0. 

PSHELL*  [16  spaces]  3 

2  .01 

2+PA 

2 

*PA  2 

1. 

2  .8333333 

0. 

PGAP 

2  0. 

0. 

400000.  0.0000  150000. 

.35 

35 

SPC 

1  291 

3456 

0. 

SPC 

1  292 

3456 

0. 

SPC 

1  293 

3456 

0. 

SPC 

1  294 

3456 

0. 

[single  point  constraints  continued] 


MPC 

5 

525 

2 

•1.00 

533 

2 

1.00 

MPC 

5 

526 

2 

•1.00 

533 

2 

1.00 

MPC 

5 

527 

2 

-1.00 

533 

2 

1.00 

MPC 

5 

528 

2 

•1.00 

533 

2 

1.00 

[multipoint  constraints  continued] 


TEMPO 

1 

970. 

TEMPD 

10 

23. 

FORCE 

2 

534 

0 

1. 

.10486 

{  . 

0. 

FORCE 

3 

534 

0 

1.50 

.10486 

v.  . 

0. 

FORCE 

4 

534 

0 

1.75 

.10486 

0. 

0. 

FORCE 

2 

534 

0 

2.00 

.10486 

0. 

0. 

FORCE 

2 

534 

0 

2.25 

.10486 

0. 

0. 

FORCE 

ENDDATA 

2 

534 

0 

2.50 

.10486 

0. 

0. 
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Appendix  B:  Mathematica  Program  for  a  Simplified 
One-Dimensional  Micromechanics  Model 
With  an  Interphase  Zone 


This  Mathematica  notebook  solves  the  constitutive  equations 
for  a  nonlinear  1 -dimensional  transverse  micromechanics  model 
containing  an  elastic -plastic  interphase  zone  between  fiber 
and  matrix. 

Input  Data  and  Clear  Previous  Run. 

Model  Parameters . 
vf  =  .35; 
fsiz  -  .07; 
zonesiz  =  .05; 
modF  =  400000.; 
modi  =  100000.; 
modM  =  100000.; 
aF  =  5.  10" (-6) ; 
al  =  10.  10"  (-6)  ; 
aM  =  10.  10" (-6)  ; 
t  =  25.  ; 
tp  =  1000  .  ; 

slys  =  125  .  ; 
hp  =  0 . ; 

Clear  [pi ,p2, sa] 


Load  arrav. 


load 

=  [0, 

10, 

20, 

30,  40 

,  50, 

60, 

70, 

80,  90 

1,  100, 

110, 

120, 

130, 

140, 

150, 

160, 

170, 

180, 

190, 

200, 

-  , 

220, 

230, 

240, 

250, 

260, 

270, 

275, 

270, 

260, 

250, 

240, 

230, 

220, 

210, 

200, 

190, 

180, 

170, 

160, 

150, 

140, 

130, 

120, 

110, 

100, 

90, 

80,  70,  60, 

50, 

40, 

30,  2C 

»,  10, 

0}  ; 

nload  =  Length [load] ; 


Define  Constitutive  Equations. 

Calculate  the  constants  a,  b,  and  c  used  in  the 
stress -strain  and  equilibrium  equations. 
a  =  fsiz; 
b  =  zonesiz  a; 

c  =  1/2  (-(2  a  +  2  b)  +  Sqrt[(2  a  +  2  b)"2  \ 

-  4  (a"2  (1  -  1/vf)  +  2  a  b  +  b"2)]); 

Stress- Strain  relations 

eF  =  uF/a  ==  sF/modF  +  aF  (t  -  tp) ; 

ell  =  (ul  -uF)/b  ==  sll/modl  +  al  (t-tp)  +  epl; 

eMl  =  (uM  -  ul)/c  ==  sMl/modM  +  aM  (t-tp); 

eI2  =  uF/a  ==  sI2/modI  +  al  (t-tp)  +  ep2; 
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eM2  =  uM/ (a+b+c)  ==  sM2/modM  +  aM  (t-tp) ; 


ql  =  sF  a  +  sI2  b  ==  sll  (a+b) ; 
q2  =  sll  ==  sMl; 

q3  =  sMl  (a+b)  +  sM2  c  ==  (a+b+c)  sa; 

Define  equation  set  and  independent  variables. 

Also,  enter .strain .hardening  equations. 
eqnset  =  {eF,eIl,eMl,eI2,eM2,ql,q2,q3} ; 
vars  =  [sF,  sll , sMl , sI2, sM2, uF,uI , uM, epl , ep2] ; 
sllys  =  slys  +  hp  Abs[pl]; 
sI2ys  =  slys  +  hp  Abs [p2] ; 

Solve  Constitutive  Equations  For  All  Possible  Cases. 

Case  1.  I sT 3  1  <  sllys  and  I  sI2 I  <  sI2vs 

solsetl  =  Solve [eqnset/ . {epi - >pl , ep2 - >p2 } , vars] ; 

fl= {solsetl [ [1,1,2] ] , solsetl  [  [1,2,2] ] , 

solsetl [[1,3,2]], solsetl [ [1,4,2] ] , solsetl  [[1,5,2]] , 

solsetl [[1,6,2]], solsetl [[1,7,2]] , solsetl [ [1,8,2] ] ,pl ,p2] ; 

Case  2.  I sll I  >  sllvs  and  I  sI2 I  <  sI2vs 

Subcase  2a.  sll  >  sllvs  and  isI2l  <  s!2ys 
solset2a  =  Solve [eqnset/ . [sll- >sllys,ep2- >p2] , vars] ; 
f 2a= [solset2a [[1,1,2]] , sllys, 

solset2a [[1,2,2]] , solset2a [[1,3,2]], sol set 2a [[1,4,2]] , 
solset2a [[1,5,2]], solset2a [[1,6,2]], solset2a  [[1,7,2]], 
solset2a [[1,8,2]] ,p2] ; 

Subcase  2b.  sll  <  -sllvs  and  lsI2l  <  sI2vs 
solset2b  =  Solve [eqnset/ . [sll ->- sllys, ep2 - >p2] , vars] ; 
f 2b= [solset2b [[1,1,2]],- sllys, 

solset2b [[1,2,2]] , sol set 2b [[1,3,2]], solset2b  [[1,4,2]] , 
solset2b [[1,5,2]] , solset2b [[1,6,2]] , solset2b [[1,7,2]], 
solset2b[ [1,8,2] ] ,p2] ; 

Case  3.  I sll I  <  sllys  and  I sI2 I  >  sI2vs 

Subcase  3a.  I sll I  <  sllys  and  sI2  >  sI2ys 

solset3a  =  Solve  [eqnset/ . [epl - >pl , sI2- >sI2ys] , vars] ; 

i 3a= [solset3a [[1,1,2]] , sol set 3a [[1,2,2]] , sol set 3a [[1,3,2]]  , 

sI2ys, solset3a [ [1,4,2] ] , sol set 3a [[1,5,2]], solset3a  [[1,6,2]] , 

solset3a[ [1,7,2]] ,pl , solset3a  [  [1 , 8, 2] ] } ; 

Subcase  3b.  [sll J  <  sllys  and  sI2  <  -sI2vs 

solset3b  =  Solve [eqnset/ . [epl- >pl , sI2- >- sI2ys] , vars] ; 

f 3b=  {solset3b [[1,1,2]] , solset3b[ [1,2,2]] , solset3b [[1,3,2]]  , 

- sI2ys, solset3b [[1,4,2]], solset3b [[1,5,2]] , solset3b [[1,6,2]] , 
solset3b [[1,7,2]] ,pl , solset3b [[1,8,2]]}; 
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Case  4.  I  sill  >  sllvs  and  lsI2  1  >  sI2ys 


Subcase  4a.  sTl  >  sllvs  and  sI2  >  sI2ys 
solset.4a  =  Solve  [eqnset/ .  {sll  -  >sllys,  sI2  -  >sI2ys}  ,vars]  ; 
f 4a= {sol set 4a [[1,1,2]] , sllys , solset4a [[1,2,2]] , sI2ys, 
solset4a [ [1,3,2]] , solset4a [[1,4,2]], solset4a [[1,5,2]] , 
solset4a [[1,6,2]], solset4a [[1,7,2]], solset4a [[1,8,2]]}; 

Subcase  4b,  sll  <  -sllvs  and  sI2  >  s!2ys 
solset4b  =  Solve [eqnset/ . {sll ->- sllys, sI2 - >sI2ys] ,vars] ; 
f 4b= {solset4b [[1,1,2]],- sllys, sol set 4b [[1,2,2]] , sI2ys, 
solset4b [[1,3,2]] , solset4b [[1,4,2]], solset4b [[1,5,2]] , 
solset 4b [ [1,6,2] ] , solset4b [ [1, 7, 2] ] , solset4b [ [1 , 8, 2] ] } ; 

Subcase  4c.  sll  >  sllvs  and  sI2  <  -sl'2ys 
solset4c  =  Solve [eqnset/ . {sll - >sllys, sI2- > -sI2ys] ,vars] ; 
f  4c= {solset 4c [[1,1,2]] , sllys , sol set 4c [[1,2,2]],- sI2ys, 
solset4c [[1,3,2]] , solset4c [[1,4,2]], solset4c [[1,5,2]] , 
solset4c [[1,6,2]] , solset4c [ [1,7,2] ] , solset4c [[1,8,2]]} ; 

Subcase  4d.  sll  <  -sllvs  and  sI2  <  -sI2vs 
solset4d  =  Solve [eqnset/ . {sll- > -sllys, sI2- >- sI2ys} ,vars] ; 
f 4d= { solset 4d [[1,1,2]], -sllys, sol set 4d [[1,2,2]], -sI2ys, 
solset 4d [ [1,3,2]] , solset 4d [[1,4,2]], solset 4d[ [1,5,2]]  , 
sol set 4d [[1,6,2]] , sol set 4d [[1,7,2]], sol set 4d [[1,8,2]]} ; 

Apply  Load  History  and  Plot  Results. 

Determine  solution  arrays.  Two  solution  arrays  are  computed, 
alvar  and  totals.  The  alvar  array  contains  the  complete 
solution  set  of  all  variables  at  each  load  step.  A 
single  element  of  alvar  is  of  the  form:  { {sF,  sll,  sMl,  sI2, 
sM2,  uF,  ul,  uM,  epl,  ep2] ,  sA}  .  Totdis  contains  only  the 
overall  displacement  of  the  structure.  A  single 
element  of  totdis  contains  {uM,  sA] . 

Pl=C . ; 
p2=0 .  ; 

n=l; 

While [n  <=  nload,  sa  =  load [ [n] ] ; 
vp  =  {0  .  ,  0  .  ,  0  .  ,  0  .  ,  0  .  ,  0  .  ,  0  .  ,  0  . ,  0  .  ,  0  . }  ; 
v  =  f  1 ; 

While  [Chop [v-vp] !  =  {0, 0, 0, 0, 0 , 0, 0, 0 , 0, 0}  , 
vp  =  v; 

v  =  Which [Abs [vp [ [2] ] ] <=sllys  &&  Abs [vp [ [4] ] ] <=sI2ys, f 1 , 

VP [ [ 2 j ]  >  sllys  &&  Abs[vp[[4]]]  <=  sI2ys,  f2a, 
vp[[2]]  <  -sllys  &&  Abs[vp[[4]]]  <=  sI2ys,  f2b, 
Abs[vp[[2]]]  <=  sllys  &&  vp[[4]]  >  sI2ys,  f3a, 

Abs  [vp  [  [2]  ]  ]  <=  sllys  &&  vp[[4]j  <  -sI2ys,  f3b, 
vp[[2]]  >  sllys  &&  VP[(4]]  >  sI2ys,  f4a, 
vp  [  [ 2 ] ]  <  -sllys  &&  vp [ [4] ]  >  sI2ys,  f4b, 
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vp[[2]]  >  sllys  &&  vp[[4]]  <  -sI2ys,  f4c, 

VP [ [2] ]  <  -sllys  &&  vp [ [ 4 ] 3  <  - sI2ys,  f4d}  ; 
pi  =  v [ [ 9 3  3  ;  P2  =  v [ [103  3  3  ; 
alvar [n] = [v, sa} ; 
totdis  [nj  =  {v  [  [83  3  ,  saj  ; 
na  =  n+1 ;  n=na3  ; 

Modify  the  total  displacement  array,  totdis.  to  contain 
ayeraq_e_s.tr.ain  after  cooldown  instead  of  total  displacement. 
Create  new  array,  es,  that  contains  the  total  average  strain 

and . applied  stress. 

length  =  a+b+c; 

es  =  Table [{ (totdis [i] [ [1] 3  -  totdis [13  [ [1] ]) /length, 
totdis [i3  [ [23 3 3 /  [i,nload}3 

Plot  the  Result. 

ListPlot [es ,  Plot Joined  - >True3 
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Appendix  C:  Mathematica  Program  for  a  Simplified 
One-Dimensional  Micromechanics  Model  With 
a  Combined  Interphase  Zone  /  Weak  Bond 


This  Mathematica  notebook  solves  the  constitutive  equations 
for  a  nonlinear  1 -dimensional  transverse  micromechanics  model 
which  models  a  portion  of  the  fiber/matrix  interface  as  if  a 
perfectly  weak  bond  exists  and  the  remainder  as  if  an 
elastic -plastic  interphase  zone  exists  between  fiber  and 
matrix. 

Input  Data  and  Clear  Previous  Run. 

Model  Parameters . 
doasiz  =  1.0 
vf  =  .35; 
fsiz  =  .07; 
zonesiz  =  .05; 
modF  =  400000.; 
modi  =  100000.; 
modM  =  100000  . ; 
aF  =  5.  10" (-6)  ; 
al  =  10.  10" (-6) ; 
aM  =  10.  10 " ( - 6)  ; 
t  =  25.  ; 
tp  =  1000 .  ; 
slys  =  125.; 
hp  =  0 .  ; 

Clear (pi, p2,sa] 

Lo-a.d_arr.ay .. 


load 

=  {0, 

10, 

20, 

30,  40 

50, 

,  60, 

70, 

80,  90 

',  100 

110, 

120, 

130, 

140, 

150, 

160, 

170, 

180, 

190, 

200, 

210, 

220, 

230, 

240, 

250, 

260, 

270, 

275, 

270, 

260, 

250, 

240, 

230, 

220, 

210, 

200, 

190, 

180, 

170, 

160, 

150, 

140, 

130, 

120, 

110, 

100, 

90, 

80,  70,  60, 

50, 

40, 

30,  20 

io 

,  0); 

nload  =  Length [load] ; 


Define  Constitutive  Equations. 

Calculate  the  constants  a.  b.  c.  and  d  used  in  the 
stress- strain  and  equilibrium  equations. 
a  =  fsiz; 
b  =  zonesiz  a; 

c  =  1/2  ( - (2  a  +  2  b'  +  Sqrt [ (2  a  +  2  b) "2  \ 

-  4  (a"2  (1  -  1/vf)  +  2  a  b  +  b"2) ] ) ; 
d  =  doasiz  a; 


121 


Stress-Strain  relations 

eF  =  uF/a  ==  sF/modF  +  aF  (t  -  tp) ; 

ell  =  (ul  -uF) /b  ==  sll/modl  +  al  (t-tp)  +  epl; 

eMl  =  (uM  -  uG) / (b+c)  ==  sMl/modM  +  aM  (t-tp); 

eI2  =  uF/a  ==  sI2/modI  +  al  (t-tp)  +  ep2; 

eM2  =  (uM  -  ul)/c  ==  sM2/modM  +  aM  (t-tp); 

eM3  =  uM/(a+b+c)  ==  sM3/modM  +  aM  (t-tp); 

Equilibrium  Equations 

ql  =  sF  a  +  sI2  b  ==  sMl  d  +  sll  (a+b-d) ; 
q2  =  sll  ==  sM2; 

q3  =  sMl  d  +  sM2  (a+b-d)  +  sM3  c  ==  (a+b+c)  sa; 

Define  equation  set  and  independent  variables. 

Also,  enter  strain  hardening  equations. 

eqnset  =  {eF,eIl,eMl,eI2,eM2,eM3,ql,q2,q3} ; 

vars  =  [sF, sll , sMl , sI2 , sM2 , sM3 , uF, ul , uM, uG, epl , ep2] ; 

sllys  =  slys  +  hp  Abs[pl] ; 

sI2ys  =  slys  +  hp  Abs[p2]; 

Solve  Constitutive  Equations  For  All  Possible  Cases. 

Case  1.  I sll I  <  sllvs  and  I s 12  I  <  sI2vs 

«  Situation  1_1 .  sMl  <  0 .  or  uG  -  uF  <  0 . 
solsetlsl  =  Solve [eqnset/ . {uG- >uF, 
epl- >pi , ep2- >p2j  , vars]  ; 

f lsl=  (solsetlsl [[1,1,2]] , solsetlsl [ [1,2,2]  ] , 
solsetlsl [[1,3,2]], solsetlsl [[1,4,2]] , solsetlsl [[1,5,2]] , 
solsetlsl [[1,6, 2]] , solsetlsl  [[1,7,2]] , 
solsetlsl [[1,8,2]], solsetlsl [[1,9,2]] , 
solsetlsl [[1,7, 2]] , pi , p2]  ; 

♦  Situation  12 .  sMl  )  0.  or  uG  ■  uF  >  0 . 
solsetls2  =  Solve [eqnset/ . [sMl - >0  .  , 
epl • : pi , ep2- >p2] , vars] ; 

f ls2= [solsetls2 [[1,1,2]] , solsetls2 [[1,2,2]],0., 
solsetls2 [[1,3,2]]  , solsetls2 [[1,4,2]] , solsetls2  [[1,5,2]] , 
sol setls2 [[1,6,2]], solsetls2  [[1,7,2]], 
solsetls2 [ [1,8,2]] , solsetls2 [[1,9,2]] ,pl ,p2] ; 


Case  2.  I sll I  >  sllvs  and  I sI2 1  <  sI2vs 


Subcase  2a.  sll  > 

sllvs  and 

1 sI2 1  < 

sI2vs 

•  Situation  2a  1 . 

sMl  <  0. 

or  uG  - 

uF  <  0. 

solset2asl  =  Solve [eqnset/ . (uG- >uF, 
sll - >sllys, ep2- >p2]  , vars] ; 
f 2asl= [solset2asl [[1,1,2]], sllys, 
solset2asl ([1,2,2]] , solset2asl [[1,3,2]] , 
solset2asl [[1,4,2]] , solset2asl ([1,5,2]]  , 
solset2asl [[1,6,2]] , solset2asl [ [1,7,2]  ]  , 
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solset2asl [[1,8,2]], solset2asl  [[1,6,2]] , 
solset2asl [ [1, 9,2] ] ,p2] ; 

»  Situation  2a_2 .  sMl  >  0.  or  uG  -  uF  >  0 . 
solset2as2  =  Solve [eqnset/ . [sMl - >0 . , 
sll - >sllys, ep2 - >p2] , vars] ; 
f 2as2=  [solset2as2 [[1,1,2]] , sllys, 0 . , 
solset2as2  [[1,2,2]],  so'lset2as2  [[1,3,2]]', 

•  solset2as2 [[1,4,2]], solset2as2 [[1,5,2]] , 
solset2as2 [[1,6,2]], solset2as2 [ [1,7,2] ] , 

'  solset2as2 [[1,8,2]], solset2as2 [[1,9,2]] ,p2] ; 

Subcase... 2b.  all,.  <_. -_sllvs_and  I  sI2 1  <  sI2vs 

•  ♦  Situation  2b_l . _ sMl  <  0  .  or  ..uG  : — uF  <_CL. 

solset2bsl  =  Solve [eqnset/ . [uG- ^uF, 

sll  sllys,  ep2 -  >p2]  ,  vars]  ; 
f 2bsl=  fsolset2bsl  [[1,1,2]],- sllys, 

£  solset2bsl [[1,2,2]] , solset2bsl [[1,3,2]] , 

solset2bsl [ [1,4,2] ] ,solset2bsl [ [1,5,2] ] , 

•  solset2bsl [[1,6,2]] , solset2bsl [ [1,7,2] ] , 

£  sclset2bsl [[1,8,2]] , solset2bsl [ [1, 6,2]  ]  , 

I  solset2bsi [ [1 , 9, 2] ] ,p2] ; 

»  Situation  2b_2 .  sMl  >  0.  or  uG  -  uF  >  0 
solset2bs2  =  Solve  [eqnset/ .  [sMl - >0 . , 

>  sll  •- sllys, ep2- >p2] , vars] ; 

\  f 2bs2= [so2  set2bs2 ([1,1,2]] , - sllys, 0., 

sc Iset 2b s 2 t  [1, 2, 2j ] , sol set 2bs2 [[1,3,2]], 
solsec2bs2  [[1,4,2]]  ,  solset2bs2  [[1,5,2]]  , 
f  solset2bs2 [[1,6,2]] , solset2bs2 [ [1 , 7 , 2] ] , 

;  solset2bs2 [[1,8,2]], solset2bs2 [[1,9,2]] , p2) ; 


Case  3.  Islll  <  silvs  and  I sI2 I  >  sI2vs 


Subcase  3a.  Islll 

<  silvs  and  sT2  > 

sI2ys 

•  Situation  3a  1 . 

sMl  <  0 .  or  uG  - 

uF  <  0 . 

»  solsetSasl  =  Solve [eqnset/ . (uG- >uF, 

epl - >pl , sI2 - >sI2ys] , vars]  ; 

l  f 3asl=  [solset3asl [[1,1,2]], solset3asl [[1,2,2]], 

solset3asl [[1,3,2]], sI2ys , sol set 3a si [[1,4,2]]  , 
solset3asl [[1,5,2]] , solset3asl [[1,6,2]]  , 
solset3asl [[1,7,2]] , solset3asl [[1,8,2]]  , 
i  solset3asl [ [1,6,2] ] ,pl , solset3asl [ [1,9,2]  ]  ]  ; 

j  •  Situation  3a_2 .  sMi  >  0.  or  uG  •  uF  >  0  ■ 

[  solset3as2  =  Solve (eqnset/ . [sMl - >0 . , 

'  epl- ;pl, sI2- >sI2ys] , vars] ; 

f  3as2=  (solset3as2  [[1,1,2]],  solset3as2  [[1,2,2]]  ,0., 
1  solset3as2 [[1,3,2]], sI2ys, solset3as2 [[1,4,2]], 

;  solset3as2 [[1,5,2]], solset3as2 [[1,6,2]], 


123 


solset3bsl  =  Solve [eqnset/ . {uG- >uF, 
epl - >pl , sI2 - > - sI2ys} , vars] ; 

f 3bsl= [solset3bsl [  [1,1,2] ] , solset3bsl [  [1,2,2] ] , 
solset3bsl [[1,3,2]], -sI2ys, solset3bsl [[1,4,2]] , 
solset3bsl [ [1,5,2] ] ,solset3bsl [ [1,6,2] ] , 
solset3bsl [[1,7,2]] , solset3bsl  [  [1 , 8, 2] ] , 
solset3bsl [ [1 , 6 , 2] ] ,pl , solset3bsl [[1,9,2]]] ; 


Situation  3b  2.  sMl  >  0,  or  uG  •  uF  >  0 


solset3bs2  =  Solve [eqnset/ . [sMl - >0 . , 
epl->pl,sI2->-sI2ys] ,vars] ; 

f 3bs2= [solset3bs2 [[1,1,2]], solset3bs2 [[1,2,2]] 
solset3bs2 [[1,3,2]],- sI2ys, solset3bs2  [[1,4,2]] 
solset3bs2 [[1,5,2]] , solset3bs2  [[1,6,2]] , 
solset3bs2 [[1,7,2]] , solset3bs2  [[1,8,2]] , 
pi, solset3bs2  [[1,9,2]]]  ; 


0. 


Case  4.  I  sill  >  sllys  and  lsI2l  >  sI2ys 


Subcase  4a.  sll  >  sllys  and  sI2  >  sI2vs 


Situation  4a  1.  sMl  <  0.  or  uG  -  uF  <  0 . 


solset4asl=Solve [eqnset/ . [uG- >uF, 
sll - >sllys, sI2 - >sI2ys] ,vars] ; 

f 4asl= [solset4asl [[1,1,2]] , sllys, solset4asl  [[1,2,2]] 
sI2ys, solset4asl [[1,3,2]], solset4asl  [[1,4,2]] , 
solset4asl [[1,5,2]], solset4asl  [[1,6,2]] , 
solset4asl [[1,7,2]], solset4asl  [[1,5,2]] , 
solset4asl [[1,8,2]], solset4asl  [[1,9,2]]]  ; 


Situation  4a  2.  sMl  >  0,  or  uG  ■  uF  >  0 . 


solset4as2=Solve  [eqnset/ .  [sMl - >0 . , 
sll- >sllys, sI2- >sI2ys] ,vars] ; 

f 4as2= [solset4as2 [[1,1,2]] , sllys, 0 . , solset4as2  [[1,2,2]] 
sI2ys, solset4as2  [[1,3,2]], solset4as2  [[1,4,2]], 
solset4as2  [[1,5,2]],  solset4as2  ['[1,6,2]]  , 
solset4as2 ([1,7,2]] , solset4as2  [[1,8,2]] , 
solset4as2 [[1,9,2]]] ; 


Subcase  4b.  sll  <  -sllvs  and  sI2  >  sI2vs 


♦  Situation  4b_l .  sMl  <  0.  or  uG  ■  uF  <  0 . 
solset4bsl=Solve [eqnset/ . [uG- >uF, 
sll- > -sllys, sI2- >sI2ys] ,vars] ; 

f 4bsl= [solset4bsl [[1,1,2]],- sllys, solset4bsl [[1,2,2]], 
sI2ys, solset4bsl [[1,3,2]], solset4bsl [[1,4,2]], 


solset4bsl [[1,5,2]] , solset4bsl [[1,6,2]]  , 
solset4bsl [ [1,7,2] ] ,solset4bsl [[1,5,2]]  , 
solset4bsl [ [1,8,2] ] , solset4bsl [[1,9,2]]}  ; 


solset4bs2=Solve [eqnset/ . [sMl->0.  , 
sll - > - sllys , sI2- >sI2ys] , vars] ; 

f 4bs2= [solset4bs2 [[1,1,2]],- sllys, 0 . , solset4bs2 [[1,2,2]] 
sI2ys, solset4bs2 [ [1,3,2] ] , solset4bs2 [[1,4,2]]  , 
solset4bs2 [ [1,5,2] ] , solset4bs2 [[1,6,2]]  , 
solset4bs2 [ [1,7,2] ] ,solset4bs2 [ [1,8,2] ] , 
solset4bs2 [[1,9,2]]}  ; 


Subcase  4c 


sllys  and  sI2 


■  sI2vs 


Situation  4c  1.  sMl  <  0 .  or  uG  -  uF  <  0, 


solset4cs2  =  Solve [eqnset/ . [sMl - >0 . , 
sll - > sllys, sI2- > - s!2ys] , vars]  ; 

f 4cs2= {soiset4cs2 [[1,1,2]] , sllys ,0., solset4cs2 [[1,2,2]], 
•  sI2ys, solset4cs2 [[1,3,2]] , solset4cs2 [[1,4,2]]  , 
selset4cs2[ [1,5,2]] , solset4cs2 [ [1, 6, 2]  ]  , 
sols«t4cs2 [ [1,7,2]] ,solset4cs2[ [1,8,2]]  , 
solset4cs2 [[1,9,2]]}  ; 


Subcase  4d 


-sllvs  and  sI2 


-sI2vs 


Situation  4d  1,  sMl  <  0 ,  or  uG _ -  uF  <  0 . 


solset4dsl  =  Solve [eqnset/ . [uG- >uF, 
sir-  *- sllys, sI2 ->- s!2ys] ,vars] ; 

f 4dsl= [solset4dsl [ [1,1,2] ] , - sllys, solset4dsl [[1,2,2]] 
-  sI2vs , solsetidsl [[1,3,2]], solset4dsl [[1,4,2]] , 
solset4dsl [[1,5,2]], solset4dsl [[1,6,2]]  , 
solset4dsl [[1,7,2]] , solset4dsl [ [1 , 5, 2] ] , 
solset4dsl [[1,8,2]], solset4dsl [[1,9,2]]}  ; 


solset4ds2  =  Solve [eqnset/ . [sMl - >0 . , 
sll • ^ -sllys, sl2 - > -sl2ysj , vars] ; 

f 4ds2= [solset4ds2 [[1,1,2]] , -sllys, 0 . , solset4ds2 [[1,2,2]] , 
- sI2ys, soiset4ds2 [ [1,3,2] ] , solset4ds2 [[1,4,2]] , 
solset4ds2 [[1,5,2]] ,solset4ds2 [ [1,6,2] ] , 
solset4ds2  [  [1,7,2]] , solset4ds2  [  [1,8,2] j , 
solset4ds2[  [1,9,2]]} ; 
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Apply  Load  History  and  Plot  Results. 


Determine  solution  arrays.  Two  solution  arrays  are  computed, 
alvar  and  totdis.  The  alvar  array  contains  the  complete 
solution  set  of  all  variables  at  each  load  step.  A  single 
element  of  alvar  is  of  the  form:  { {sF,  sll,  sMl,  sI2,  sM2, 
uF,  ul,  uM,  epl,  ep2}  ,  sA] .  Totdis  contains  only  the  overall 
displacement  of  the  structure.  A  single  element  of  totdis 
contains  [uM,  sA} . 


pl=0 . ; 
p2=0 . ; 

n=l ; 

While [n  <=  nload,  sa  =  load[[n]]; 
vp  =  [0.,0.,0.,0.,0.,0.f0.,0.,0.,0.,0.,0.}; 
v  =  fls2; 

While [Chop [v-vp]  !  =  [0,0,0,0,0,0,0, 0,0, 0,0,0}  , 
vp  =  v; 

uGF  =  vp[  [10]  ]  -vp[  [7]  ]  ; 

v  =  Which [Abs [vp [ [2] ] ] <=sllys  &&  Abs [vp [ [4] ] ] <=sI2ys 
&&  (vp [ [3] ] <0 .  ||  uGF<0) ,  flsl, 

Abs[vp[[2]]]  <=  sllys  &&  Abs[vp[[4]]]  <=  sI2ys 
&&  (vp  [  [3]  ]  >0 .  ||  uGF>0)  ,  f  ls2, 
vp  [  [2]  ]  >  sllys  &.Sc  Abs[vp[[4]]]  <=  sI2ys 
&&  (vp [ [3] ] <0 .  ||  uGF<0) , f 2asl , 

VP  [  [2] ]  >  sllys  &&  Abs [vp [ [4] ] ]  <=  sI2ys 

&&  (vp [ [3] ] >0 .  ||  uGF>0) , f 2as2 , 

vp [  [2] ]  <  -sllys  &&  Abs[vp[[4]j]  <=  sI2ys 

&&  (vp  [  [3] ] <0 .  ||  uGF < 0)  ,  f 2bsl , 

vp [  [2] ]  <  -sllys  &&  Abs[vp[[4]]]  <=  sI2ys 

&&  (vp [ [3] ] >0 .  ||  uGF>0) ,  f 2bs2 , 

Abs(vp[[2]]]  <=  sllys  &&  vp  [  [4] ]  >  sI2ys 
&&  (vp [ [3] ] <0 .  ||  uGF<0) ,  f 3asl , 

Abs[vp[[2]]]  <=  sllys  &&  vp[[4]j  >  sI2ys 
U&  (vp  [  [3] ] >0 .  ||  uGF>0) ,  f 3as2 , 

Abs[vp[[2]j]  <=  sllys  &&  vp  [  [ 4 ] ]  <  -sI2ys 
&&  (vp [  [3] ] <0 .  ||  uGF<0) ,  f 3bsl , 

Abs[vp[[2]]]  <=  sllys  &&  vp[ [4] ]  <  -sI2ys 
&&  (vp [ [3] ] >0 .  ||  uGF>0) ,  f 3bs2 , 
vp [ [ 2 ] ]  >  sllys  &&  vp[[4]]  >  sI2ys 
&&  (vp [ [3] ] <0 .  ||  uGF <0) ,  f 4asl , 
vp [  [2] ]  >  sllys  &&  vp[[4J]  >  sI2ys 
&&  (vp  [  [3] ] >0 .  ||  uGF>0) ,  f  4as2 , 
vp [ [2] ]  <  -sllys  &&  vp[[4]]  >  sI2ys 
&&  (vp  [  [3] ] <0 .  ||  uGF<0) ,  f 4bsl , 

VP[[2]]  <  -sllys  &&  vp[[4]]  >  sI2ys 
&&  (vp [ [3] ] >0 .  II  uGF>0) ,  f 4bs2 , 
vp [  [2] ]  >  sllys  &&  vp  [  [4] ]  <  -sI2ys 
&&  (vp [ [3] ] <0 .  ||  uGF<0) ,  f 4csl , 
vp[[2]]  >  sllys  &&  vp[[4]]  <  -sI2ys 
&&  (vp [ [3] ] >0 .  ||  uGF>0) ,  f 4cs2, 
vp  [  [2] ]  <  -sllys  &&  vp[[4]]  <  -sI2ys 


126 


&&  (vp [ [3] ] <0 .  II  uGF<0) ,  f 4dsl , 
vp[[2]]  <  -sllys  &&  vp[[4]]  <  -sI2ys 
&&  (vp [ [3] ] >0 .  ||  uGF>0) ,  f4ds2] ; 
pi  =  v[  [11]  ]  ;  p2  =  v  [  [12]  ]  ]  ; 
alvar [n] = {v, sa] ; 
totdis [n] = {v [ [9] ] , sa] ; 
na  =  n+1 ;  n=na] ; 

MQdif.v_the  total  displacement  array,  totdis.  to  contain 
average  strain  after  cooldown  instead  of  total  displacement. 
Crfiat.e_new.„arrav._eJaJ_that.  contains  the  total,  a v_eraa.e„ strain 
and  applied  stress. 
length  =  a+b+c; 

es  =  Table [{ (totdis [i]  [[1]]  -  totdis  (1]  [  [1] ]) /length, 
totdis [i] [ [2] ]} ,  {i,nload}] 

Bl.Q.t-.the...Result . 

ListPlot[es,  PlotJoined- >True] 
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